MODULE INIT

! Allocate numerous arrays and perform miscellaneous initializations

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE MESH_POINTERS
USE GLOBAL_CONSTANTS
USE OUTPUT_DATA
USE TRAN
USE MEMORY_FUNCTIONS, ONLY : CHKMEMERR
USE DEVICE_VARIABLES

IMPLICIT NONE

PRIVATE

PUBLIC INITIALIZE_MESH_VARIABLES_1,INITIALIZE_MESH_VARIABLES_2,INITIALIZE_GLOBAL_VARIABLES, &
       OPEN_AND_CLOSE,INITIAL_NOISE,UVW_INIT,INITIALIZE_DEVICES,INITIALIZE_PROFILES,REASSIGN_WALL_CELLS

CONTAINS


SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)

USE PHYSICAL_FUNCTIONS, ONLY: GET_VISCOSITY,GET_SPECIFIC_GAS_CONSTANT,GET_SPECIFIC_HEAT,LES_FILTER_WIDTH_FUNCTION
USE MEMORY_FUNCTIONS, ONLY: COMPUTE_ONE_D_STORAGE_DIMENSIONS,COMPUTE_PARTICLE_STORAGE_DIMENSIONS,COMPUTE_WALL_STORAGE_DIMENSIONS,&
COMPUTE_CFACE_STORAGE_DIMENSIONS
USE GEOMETRY_FUNCTIONS, ONLY: ASSIGN_HT3D_WALL_INDICES
USE RADCONS, ONLY: UIIDIM
USE CONTROL_VARIABLES
USE MATH_FUNCTIONS, ONLY:EVALUATE_RAMP
INTEGER :: N,NN,I,J,K,IW,IC,SURF_INDEX,IOR,IERR,IZERO,II,JJ,KK,I_OBST
REAL(EB), INTENT(IN) :: DT
INTEGER, INTENT(IN) :: NM
REAL(EB) :: MU_N,ZZ_GET(1:N_TRACKED_SPECIES),CS,DELTA,DUMMY=0._EB,INTEGRAL,TEMP,ZSW
INTEGER, POINTER :: IBP1, JBP1, KBP1,IBAR, JBAR, KBAR, N_EDGES
INTEGER, TARGET  :: NULL_TARGET=NULL_BOUNDARY,DEFAULT_SURF_INDEX_TARGET=0
REAL(EB),POINTER :: XS,XF,YS,YF,ZS,ZF
TYPE (INITIALIZATION_TYPE), POINTER :: IN=>NULL()
TYPE (VENTS_TYPE), POINTER :: VT=>NULL()
TYPE (OBSTRUCTION_TYPE), POINTER :: OB=>NULL()
TYPE (WALL_TYPE), POINTER :: WC=>NULL()
TYPE (SURFACE_TYPE), POINTER :: SF=>NULL()
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE (RAMPS_TYPE), POINTER :: RP=>NULL()
TYPE (MULTIPLIER_TYPE), POINTER :: MR=>NULL()


IERR = 0
M => MESHES(NM)
IBP1 =>M%IBP1
JBP1 =>M%JBP1
KBP1 =>M%KBP1
IBAR =>M%IBAR
JBAR =>M%JBAR
KBAR =>M%KBAR
N_EDGES=>M%N_EDGES
XS=>M%XS
YS=>M%YS
ZS=>M%ZS
XF=>M%XF
YF=>M%YF
ZF=>M%ZF

ALLOCATE(M%EXTERNAL_WALL(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','EXTERNAL_WALL',IZERO)

ALLOCATE(M%RHO(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO',IZERO)
ALLOCATE(M%RHOS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHOS',IZERO)
M%RHOS = RHOA
ALLOCATE(M%TMP(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP',IZERO)
ALLOCATE(M%U(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','U',IZERO)
ALLOCATE(M%V(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','V',IZERO)
ALLOCATE(M%W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','W',IZERO)
ALLOCATE(M%US(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','US',IZERO)
ALLOCATE(M%VS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','VS',IZERO)
ALLOCATE(M%WS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WS',IZERO)
ALLOCATE(M%FVX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVX',IZERO)
ALLOCATE(M%FVY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVY',IZERO)
ALLOCATE(M%FVZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVZ',IZERO)
ALLOCATE(M%FVX_B(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVX_B',IZERO) ; M%FVX_B=0._EB
ALLOCATE(M%FVY_B(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVY_B',IZERO) ; M%FVY_B=0._EB
ALLOCATE(M%FVZ_B(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVZ_B',IZERO) ; M%FVZ_B=0._EB
ALLOCATE(M%H(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','H',IZERO)
ALLOCATE(M%HS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','HS',IZERO)
IF (EXTERNAL_BOUNDARY_CORRECTION) THEN
   ALLOCATE(M%H_PRIME(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','H_PRIME',IZERO)
ENDIF
ALLOCATE(M%KRES(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','KRES',IZERO)
ALLOCATE(M%DDDT(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DDDT',IZERO)
ALLOCATE(M%D(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','D',IZERO)
ALLOCATE(M%DS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DS',IZERO)
ALLOCATE(M%MU(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MU',IZERO)
ALLOCATE(M%MU_DNS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MU_DNS',IZERO)
IF (CHECK_VN) THEN
   ALLOCATE(M%D_Z_MAX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','D_Z_MAX',IZERO)
   M%D_Z_MAX=0._EB
ENDIF
ALLOCATE(M%STRAIN_RATE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','STRAIN_RATE',IZERO)
M%STRAIN_RATE = 0._EB
SELECT CASE(TURB_MODEL)
   CASE (CONSMAG,DYNSMAG)
      ALLOCATE(M%CSD2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
      CALL ChkMemErr('INIT','CS',IZERO)
END SELECT
IF (OUTPUT_CHEM_IT) THEN
   ALLOCATE(M%CHEM_SUBIT(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','CHEM_SUBIT',IZERO)
   M%CHEM_SUBIT = 0._EB
ENDIF
ALLOCATE(M%Q(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','Q',IZERO)

ALLOCATE(M%MIX_TIME(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MIX_TIME',IZERO)
M%MIX_TIME=DT

IF (POTENTIAL_TEMPERATURE_CORRECTION) THEN
   ALLOCATE(M%PR_T(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','PR_T',IZERO)
   M%PR_T=PR
ENDIF

IF (TRI_MODEL) THEN
   ALLOCATE(M%TRI_COR(0:IBP1,0:JBP1,0:KBP1,NUMBER_SPECTRAL_BANDS),STAT=IZERO)
   CALL ChkMemErr('INIT','TRI_COR',IZERO)
   M%TRI_COR=1._EB
   ALLOCATE(M%TMP_FLAME(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','TMP_FLAME',IZERO)
   M%TMP_FLAME=TMPA
ENDIF

IF (FLAME_INDEX_MODEL) THEN
   ALLOCATE(M%FLAME_INDEX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','FLAME_INDEX',IZERO)
   M%FLAME_INDEX=-1._EB
ENDIF

! Background pressure, temperature, density as a function of height (Z coordinate)

ALLOCATE(M%PRESSURE_ZONE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','PRESSURE_ZONE',IZERO) ; M%PRESSURE_ZONE = -1

ALLOCATE(M%P_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','P_0',IZERO)
ALLOCATE(M%TMP_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP_0',IZERO)
ALLOCATE(M%RHO_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO_0',IZERO)

! Allocate species arrays

IF (.NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE( M%ZZ(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ZZ',IZERO)
   M%ZZ = 0._EB
   ALLOCATE(M%ZZS(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ZZS',IZERO)
   M%ZZS = 0._EB
   ALLOCATE(M%DEL_RHO_D_DEL_Z(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DEL_RHO_D_DEL_Z',IZERO)
   M%DEL_RHO_D_DEL_Z = 0._EB
   IF (REAC_SOURCE_CHECK) THEN
      ALLOCATE( M%REAC_SOURCE_TERM(0:IBP1,0:JBP1,0:KBP1,N_TRACKED_SPECIES),STAT=IZERO)
      CALL ChkMemErr('INIT','REAC_SOURCE_TERM',IZERO)
      M%REAC_SOURCE_TERM = 0._EB
      ALLOCATE( M%Q_REAC(0:IBP1,0:JBP1,0:KBP1,N_REACTIONS),STAT=IZERO)
      CALL ChkMemErr('INIT','Q_REAC',IZERO)
      M%Q_REAC = 0._EB
   ENDIF
ENDIF

ALLOCATE(M%RSUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RSUM',IZERO)
M%RSUM = RSUM0

! Allocate scalar face values

ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','FX',IZERO)
ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','FY',IZERO)
ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','FZ',IZERO)
M%FX = 0._EB
M%FY = 0._EB
M%FZ = 0._EB

! Allocate storage for scalar total fluxes

IF (STORE_SPECIES_FLUX) THEN
   ALLOCATE( M%ADV_FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ADV_FX',IZERO)
   ALLOCATE( M%ADV_FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ADV_FY',IZERO)
   ALLOCATE( M%ADV_FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ADV_FZ',IZERO)
   ALLOCATE( M%DIF_FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FX',IZERO)
   ALLOCATE( M%DIF_FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FY',IZERO)
   ALLOCATE( M%DIF_FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FZ',IZERO)
   ALLOCATE( M%DIF_FXS(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FX',IZERO)
   ALLOCATE( M%DIF_FYS(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FY',IZERO)
   ALLOCATE( M%DIF_FZS(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FZ',IZERO)
   M%ADV_FX = 0._EB
   M%ADV_FY = 0._EB
   M%ADV_FZ = 0._EB
   M%DIF_FX = 0._EB
   M%DIF_FY = 0._EB
   M%DIF_FZ = 0._EB
   M%DIF_FXS = 0._EB
   M%DIF_FYS = 0._EB
   M%DIF_FZS = 0._EB
ENDIF

! Allocate water mass arrays if sprinklers are present

IF (N_LP_ARRAY_INDICES>0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%QR_W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR_W',IZERO)
   M%QR_W = 0._EB
   ALLOCATE(M%AVG_DROP_DEN(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_DEN',IZERO)
   M%AVG_DROP_DEN=0._EB
   ALLOCATE(M%AVG_DROP_AREA(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_AREA',IZERO)
   M%AVG_DROP_AREA=0._EB
   ALLOCATE(M%AVG_DROP_TMP(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_TMP',IZERO)
   M%AVG_DROP_TMP=TMPM
   ALLOCATE(M%AVG_DROP_RAD(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_RAD',IZERO)
   M%AVG_DROP_RAD=0._EB
ENDIF

IF (N_LP_ARRAY_INDICES>0 .OR. N_REACTIONS>0 .OR. ANY(SPECIES_MIXTURE%DEPOSITING) .OR. &
    ANY(SPECIES_MIXTURE%CONDENSATION_SMIX_INDEX>0)) THEN
   ALLOCATE(M%D_SOURCE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','D_SOURCE',IZERO)
   M%D_SOURCE = 0._EB
   ALLOCATE(M%M_DOT_PPP(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','M_DOT_PPP',IZERO)
   M%M_DOT_PPP=0._EB
ENDIF

! If radiation absorption desired allocate arrays

IF (.NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%CHI_R(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR',IZERO)
   IF (N_REACTIONS>0) THEN
      M%CHI_R = MAX(0._EB,RADIATIVE_FRACTION*EVALUATE_RAMP(T_BEGIN,0._EB,REACTION(1)%RAMP_CHI_R_INDEX))
   ELSE
      M%CHI_R = MAX(0._EB,RADIATIVE_FRACTION)
   ENDIF
   ALLOCATE(M%QR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR',IZERO)
   M%QR = 0._EB
   ALLOCATE(M%KAPPA_GAS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','KAPPA_GAS',IZERO)
   M%KAPPA_GAS = 0._EB
   ALLOCATE(M%UII(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','UII',IZERO)
   M%UII = 0._EB
ELSE
   ALLOCATE(M%QR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR',IZERO)
   M%QR = 0._EB
ENDIF

! Work arrays

ALLOCATE(M%WORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK1',IZERO)
ALLOCATE(M%WORK2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK2',IZERO)
ALLOCATE(M%WORK3(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK3',IZERO)
ALLOCATE(M%WORK4(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK4',IZERO)
ALLOCATE(M%WORK5(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK5',IZERO)
ALLOCATE(M%WORK6(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK6',IZERO)
ALLOCATE(M%WORK7(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK7',IZERO)
ALLOCATE(M%WORK8(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK8',IZERO)
ALLOCATE(M%WORK9(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK9',IZERO)

! for testing corrections to background pressure due to variable specific heat
IF (STORE_DIVERGENCE_CORRECTION) THEN
   ALLOCATE(M%DCOR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','DCOR',IZERO)
   M%DCOR=0._EB
ENDIF

ALLOCATE(M%IWORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','IWORK1',IZERO)
ALLOCATE(M%SCALAR_WORK1(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','SCALAR_WORK1',IZERO)
ALLOCATE(M%SCALAR_WORK2(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','SCALAR_WORK2',IZERO)
ALLOCATE(M%SCALAR_WORK3(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','SCALAR_WORK3',IZERO)
ALLOCATE(M%SCALAR_WORK4(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','SCALAR_WORK4',IZERO)
M%IWORK1=0
M%SCALAR_WORK1=0._EB
M%SCALAR_WORK2=0._EB
M%SCALAR_WORK3=0._EB
M%SCALAR_WORK4=0._EB

! Boundary file patch counter

ALLOCATE(M%INC(-3:3,0:M%N_OBST),STAT=IZERO)
CALL ChkMemErr('INIT','INC',IZERO)

IF (STRATIFICATION .AND. .NOT.EVACUATION_ONLY(NM)) THEN

   ! Compute the atmospheric pressure profile ramp using the specified temperature ramp

   RP => RAMPS(I_RAMP_P0_Z)
   INTEGRAL = 0._EB
   IF (HVAC_SOLVE) THEN
      ZSW = MIN(ZS_MIN-DZS_MAX,NODE_Z_MIN)
   ELSE
      ZSW = ZS_MIN
   ENDIF

   DO K=0,RP%NUMBER_INTERPOLATION_POINTS+1
      TEMP = TMPA*RAMPS(I_RAMP_TMP0_Z)%INTERPOLATED_DATA(K)
      INTEGRAL = INTEGRAL + (GVEC(3)/(RSUM0*TEMP))/RP%RDT
      RP%INTERPOLATED_DATA(K) = P_INF*EXP(GVEC(3)*(ZSW-GROUND_LEVEL)/(RSUM0*TMPA))*EXP(INTEGRAL)
   ENDDO

   ! Populate the cell-centered background temperature and pressure

   DO K=0,M%KBP1
      M%TMP_0(K) = TMPA*EVALUATE_RAMP(M%ZC(K),DUMMY,I_RAMP_TMP0_Z)
      M%P_0(K)   = EVALUATE_RAMP(M%ZC(K),DUMMY,I_RAMP_P0_Z)
   ENDDO
ELSE

   M%TMP_0(:) = TMPA
   M%P_0(:)   = P_INF

ENDIF

! Initialize density profile

DO K=0,M%KBP1
   M%RHO_0(K) = M%P_0(K)/(M%TMP_0(K)*RSUM0)
ENDDO

! Initialize various time step variables

DT_INITIAL = DT

! Initialize major arrays

DO K=0,M%KBP1
   M%RHO(:,:,K) = M%RHO_0(K)
   M%TMP(:,:,K) = M%TMP_0(K)
   M%U(:,:,K)   = U0*EVALUATE_RAMP(T_BEGIN,DUMMY,I_RAMP_U0_T)*EVALUATE_RAMP(M%ZC(K),DUMMY,I_RAMP_U0_Z)
   M%V(:,:,K)   = V0*EVALUATE_RAMP(T_BEGIN,DUMMY,I_RAMP_V0_T)*EVALUATE_RAMP(M%ZC(K),DUMMY,I_RAMP_V0_Z)
   M%W(:,:,K)   = W0*EVALUATE_RAMP(T_BEGIN,DUMMY,I_RAMP_W0_T)*EVALUATE_RAMP(M%ZC(K),DUMMY,I_RAMP_W0_Z)
ENDDO

IF (ANY(MEAN_FORCING) .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   DO K=0,M%KBP1
      DO J=0,M%JBP1
         DO I=0,M%IBAR
            IF ( .NOT.(M%MEAN_FORCING_CELL(I,J,K) .AND. M%MEAN_FORCING_CELL(I+1,J,K)) ) M%U(I,J,K)=0._EB
         ENDDO
      ENDDO
   ENDDO
   DO K=0,M%KBP1
      DO J=0,M%JBAR
         DO I=0,M%IBP1
            IF ( .NOT.(M%MEAN_FORCING_CELL(I,J,K) .AND.  M%MEAN_FORCING_CELL(I,J+1,K)) ) M%V(I,J,K)=0._EB
         ENDDO
      ENDDO
   ENDDO
   DO K=0,M%KBAR
      DO J=0,M%JBP1
         DO I=0,M%IBP1
            IF ( .NOT.(M%MEAN_FORCING_CELL(I,J,K) .AND.  M%MEAN_FORCING_CELL(I,J,K+1)) ) M%W(I,J,K)=0._EB
         ENDDO
      ENDDO
   ENDDO
ENDIF

M%US    = M%U
M%VS    = M%V
M%WS    = M%W
M%FVX   = 0._EB
M%FVY   = 0._EB
M%FVZ   = 0._EB
M%H     = H0
M%HS    = H0
M%KRES  = 0._EB
IF (EXTERNAL_BOUNDARY_CORRECTION) M%H_PRIME=0._EB

M%DDDT  = 0._EB
M%D     = 0._EB
M%DS    = 0._EB
M%Q  = 0._EB
IF (EVACUATION_ONLY(NM)) THEN
   M%U  = 0._EB
   M%V  = 0._EB
   M%W  = 0._EB
   M%US = 0._EB
   M%VS = 0._EB
   M%WS = 0._EB
   M%H  = 0._EB
   M%HS = 0._EB
ENDIF

! Viscosity

ZZ_GET(1:N_TRACKED_SPECIES) = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
CALL GET_VISCOSITY(ZZ_GET,MU_N,TMPA)
M%MU = MU_N

CS = C_SMAGORINSKY
IF (EVACUATION_ONLY(NM)) CS=0.9_EB
SELECT CASE(TURB_MODEL)
   CASE(CONSMAG,DYNSMAG)
      DO K=0,KBP1
         DO J=0,JBP1
            DO I=0,IBP1
               DELTA = LES_FILTER_WIDTH_FUNCTION(M%DX(I),M%DY(J),M%DZ(K))
               M%CSD2(I,J,K) = (CS*DELTA)**2
            ENDDO
         ENDDO
      ENDDO
END SELECT

! Initialize mass fraction arrays

IF (.NOT.EVACUATION_ONLY(NM)) THEN
   DO N=1,N_TRACKED_SPECIES
      M%ZZ(:,:,:,N)  = SPECIES_MIXTURE(N)%ZZ0
      M%ZZS(:,:,:,N) = SPECIES_MIXTURE(N)%ZZ0
   ENDDO
   DO N=N_TRACKED_SPECIES+1,N_TOTAL_SCALARS
      M%ZZ(:,:,:,N)  = INITIAL_UNMIXED_FRACTION
      M%ZZS(:,:,:,N) = INITIAL_UNMIXED_FRACTION
   ENDDO
ENDIF

! Over-ride default ambient conditions with user-prescribed INITializations

DO N=1,N_INIT
   IF (EVACUATION_ONLY(NM)) CYCLE
   IN => INITIALIZATION(N)
   DO K=0,KBP1
      DO J=0,JBP1
         DO I=0,IBP1
            IF (M%XC(I) > IN%X1 .AND. M%XC(I) < IN%X2 .AND. &
                M%YC(J) > IN%Y1 .AND. M%YC(J) < IN%Y2 .AND. &
                M%ZC(K) > IN%Z1 .AND. M%ZC(K) < IN%Z2) THEN
               M%TMP(I,J,K)            = IN%TEMPERATURE
               M%RHO(I,J,K)            = IN%DENSITY
               M%RHOS(I,J,K)           = IN%DENSITY
               M%ZZ(I,J,K,1:N_TRACKED_SPECIES) = IN%MASS_FRACTION(1:N_TRACKED_SPECIES)
               M%ZZS(I,J,K,1:N_TRACKED_SPECIES) = IN%MASS_FRACTION(1:N_TRACKED_SPECIES)
               IF (IN%ADJUST_DENSITY)     M%RHO(I,J,K) = M%RHO(I,J,K)*M%P_0(K)/P_INF
               IF (IN%ADJUST_DENSITY)     M%RHOS(I,J,K) = M%RHOS(I,J,K)*M%P_0(K)/P_INF
               IF (IN%ADJUST_TEMPERATURE) M%TMP(I,J,K) = M%TMP(I,J,K)*M%P_0(K)/P_INF
            ENDIF
         ENDDO
      ENDDO
   ENDDO
ENDDO

! Compute molecular weight term RSUM=R0*SUM(Y_i/M_i)

IF (.NOT.EVACUATION_ONLY(NM)) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            ZZ_GET(1:N_TRACKED_SPECIES) = M%ZZ(I,J,K,1:N_TRACKED_SPECIES)
            CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,M%RSUM(I,J,K))
         ENDDO
      ENDDO
   ENDDO
ENDIF

! Allocate and Initialize Mesh-Dependent Radiation Arrays

M%QR    = 0._EB
IF (.NOT.EVACUATION_ONLY(NM)) M%UII = 4._EB*SIGMA*TMPA4
M%ANGLE_INC_COUNTER = 0
M%RAD_CALL_COUNTER  = 0
IF (RADIATION .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%UIID(0:M%IBP1,0:M%JBP1,0:M%KBP1,1:UIIDIM),STAT=IZERO)
   CALL ChkMemErr('INIT','UIID',IZERO)
   M%UIID = 4._EB*SIGMA*TMPA4/REAL(UIIDIM,EB)
ENDIF

! General work arrays

M%WORK1 = 0._EB
M%WORK2 = 0._EB
M%WORK3 = 0._EB
M%WORK4 = 0._EB
M%WORK5 = 0._EB
M%WORK6 = 0._EB
IF (.NOT.EVACUATION_ONLY(NM)) M%WORK7 = 0._EB

! Allocate ONE_D data storage array and compute the dimensions of its components

ALLOCATE(M%WALL_STORAGE(0:N_SURF))
ALLOCATE(M%CFACE_STORAGE(0:N_SURF))

DO N=0,N_SURF
   CALL COMPUTE_ONE_D_STORAGE_DIMENSIONS(N)
   CALL COMPUTE_WALL_STORAGE_DIMENSIONS(N)
   CALL COMPUTE_CFACE_STORAGE_DIMENSIONS(N)
ENDDO

! Allocate lagrangian particle storage array and compute the dimensions of its components

M%NLP = 0
M%NLPDIM = 1000
M%PARTICLE_TAG = NM

IF (PARTICLE_FILE .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%LAGRANGIAN_PARTICLE(M%NLPDIM),STAT=IZERO)
   CALL ChkMemErr('INIT','PARTICLE',IZERO)
   ALLOCATE(M%PARTICLE_STORAGE(N_LAGRANGIAN_CLASSES))
   DO N=1,N_LAGRANGIAN_CLASSES
      CALL COMPUTE_PARTICLE_STORAGE_DIMENSIONS(N)
   ENDDO
ENDIF

! Allocate wall cell arrays. This array will grow as needed.

M%N_WALL_CELLS_DIM = 1000

ALLOCATE(M%WALL(0:M%N_WALL_CELLS_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','WALL',IZERO)

M%WALL(0)%BOUNDARY_TYPE => NULL_TARGET
M%WALL(0)%SURF_INDEX    => DEFAULT_SURF_INDEX_TARGET

! Allocate arrays that are dimensioned by the number of external wall cells

NOT_EVAC_IF_1: IF (.NOT.EVACUATION_ONLY(NM)) THEN

   ALLOCATE(M%D_CORR(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
   CALL ChkMemErr('INIT','D_CORR',IZERO)
   ALLOCATE(M%DS_CORR(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
   CALL ChkMemErr('INIT','DS_CORR',IZERO)
   M%D_CORR = 0._EB
   M%DS_CORR = 0._EB
   ALLOCATE(M%UVW_SAVE(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
   CALL ChkMemErr('INIT','UVW_SAVE',IZERO)
   M%UVW_SAVE = 0._EB

   ALLOCATE(M%U_GHOST(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
   CALL ChkMemErr('INIT','U_GHOST',IZERO)
   ALLOCATE(M%V_GHOST(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
   CALL ChkMemErr('INIT','V_GHOST',IZERO)
   ALLOCATE(M%W_GHOST(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
   CALL ChkMemErr('INIT','W_GHOST',IZERO)
   M%U_GHOST = 0._EB
   M%V_GHOST = 0._EB
   M%W_GHOST = 0._EB
ENDIF NOT_EVAC_IF_1

! Allocate arrays for turbulent inflow boundary conditions (experimental)

VENT_LOOP: DO N=1,M%N_VENT
   VT => M%VENTS(N)
   EDDY_IF: IF (VT%N_EDDY>0) THEN
      SELECT CASE(ABS(VT%IOR))
         CASE(1)
            ALLOCATE(VT%U_EDDY(VT%J1+1:VT%J2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','U_EDDY',IZERO)
            ALLOCATE(VT%V_EDDY(VT%J1+1:VT%J2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','V_EDDY',IZERO)
            ALLOCATE(VT%W_EDDY(VT%J1+1:VT%J2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','W_EDDY',IZERO)
         CASE(2)
            ALLOCATE(VT%U_EDDY(VT%I1+1:VT%I2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','U_EDDY',IZERO)
            ALLOCATE(VT%V_EDDY(VT%I1+1:VT%I2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','V_EDDY',IZERO)
            ALLOCATE(VT%W_EDDY(VT%I1+1:VT%I2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','W_EDDY',IZERO)
         CASE(3)
            ALLOCATE(VT%U_EDDY(VT%I1+1:VT%I2,VT%J1+1:VT%J2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','U_EDDY',IZERO)
            ALLOCATE(VT%V_EDDY(VT%I1+1:VT%I2,VT%J1+1:VT%J2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','V_EDDY',IZERO)
            ALLOCATE(VT%W_EDDY(VT%I1+1:VT%I2,VT%J1+1:VT%J2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','W_EDDY',IZERO)
      END SELECT
      ALLOCATE(VT%X_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','X_EDDY',IZERO)
      ALLOCATE(VT%Y_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','Y_EDDY',IZERO)
      ALLOCATE(VT%Z_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','Z_EDDY',IZERO)
      ALLOCATE(VT%CU_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','CU_EDDY',IZERO)
      ALLOCATE(VT%CV_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','CV_EDDY',IZERO)
      ALLOCATE(VT%CW_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','CW_EDDY',IZERO)
      VT%U_EDDY=0._EB
      VT%V_EDDY=0._EB
      VT%W_EDDY=0._EB
      VT%X_EDDY=0._EB
      VT%Y_EDDY=0._EB
      VT%Z_EDDY=0._EB
      VT%CU_EDDY=0._EB
      VT%CV_EDDY=0._EB
      VT%CW_EDDY=0._EB
   ENDIF EDDY_IF
ENDDO VENT_LOOP

! Allocate wall and edge index arrays

ALLOCATE(M%WALL_INDEX_HT3D(0:CELL_COUNT(NM),-3:3),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_INDEX_HT3D',IZERO)
M%WALL_INDEX_HT3D = 0

ALLOCATE(M%EDGE_INDEX(1:12,0:CELL_COUNT(NM)),STAT=IZERO)
CALL ChkMemErr('INIT','EDGE_INDEX',IZERO)
M%EDGE_INDEX = 0

ALLOCATE(M%U_EDGE_Y(0:CELL_COUNT(NM)),STAT=IZERO) ; CALL ChkMemErr('INIT','U_EDGE_Y',IZERO) ; M%U_EDGE_Y = -1.E6_EB
ALLOCATE(M%U_EDGE_Z(0:CELL_COUNT(NM)),STAT=IZERO) ; CALL ChkMemErr('INIT','U_EDGE_Z',IZERO) ; M%U_EDGE_Z = -1.E6_EB
ALLOCATE(M%V_EDGE_X(0:CELL_COUNT(NM)),STAT=IZERO) ; CALL ChkMemErr('INIT','V_EDGE_X',IZERO) ; M%V_EDGE_X = -1.E6_EB
ALLOCATE(M%V_EDGE_Z(0:CELL_COUNT(NM)),STAT=IZERO) ; CALL ChkMemErr('INIT','V_EDGE_Z',IZERO) ; M%V_EDGE_Z = -1.E6_EB
ALLOCATE(M%W_EDGE_X(0:CELL_COUNT(NM)),STAT=IZERO) ; CALL ChkMemErr('INIT','W_EDGE_X',IZERO) ; M%W_EDGE_X = -1.E6_EB
ALLOCATE(M%W_EDGE_Y(0:CELL_COUNT(NM)),STAT=IZERO) ; CALL ChkMemErr('INIT','W_EDGE_Y',IZERO) ; M%W_EDGE_Y = -1.E6_EB

! Set up WALL for external boundaries of the current mesh

M%N_WALL_CELLS = 0

DO IOR=1,-1,-2
   IF (IOR== 1) I = 0
   IF (IOR==-1) I = IBP1
   DO K=1,KBAR
      DO J=1,JBAR
         M%N_WALL_CELLS = M%N_WALL_CELLS + 1
         CALL INIT_WALL_CELL(NM,I,J,K,0,M%N_WALL_CELLS,IOR,DEFAULT_SURF_INDEX,IERR,T_BEGIN) ; IF (IERR>0) RETURN
      ENDDO
   ENDDO
ENDDO

DO IOR=2,-2,-4
   IF (IOR== 2) J = 0
   IF (IOR==-2) J = JBP1
   DO K=1,KBAR
      DO I=1,IBAR
         M%N_WALL_CELLS = M%N_WALL_CELLS + 1
         CALL INIT_WALL_CELL(NM,I,J,K,0,M%N_WALL_CELLS,IOR,DEFAULT_SURF_INDEX,IERR,T_BEGIN) ; IF (IERR>0) RETURN
      ENDDO
   ENDDO
ENDDO

NOT_EVAC_IF_2: IF (.NOT.EVACUATION_ONLY(NM)) THEN
DO IOR=3,-3,-6
   IF (IOR== 3) K = 0
   IF (IOR==-3) K = KBP1
   DO J=1,JBAR
      DO I=1,IBAR
         M%N_WALL_CELLS = M%N_WALL_CELLS + 1
         CALL INIT_WALL_CELL(NM,I,J,K,0,M%N_WALL_CELLS,IOR,DEFAULT_SURF_INDEX,IERR,T_BEGIN) ; IF (IERR>0) RETURN
      ENDDO
   ENDDO
ENDDO
ENDIF NOT_EVAC_IF_2

! Go through all obstructions and decide which cell faces ought to be given a wall cell index and initialized

M%N_INTERNAL_WALL_CELLS = 0

OBST_LOOP_2: DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)

   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I1+1
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (I==1) CYCLE
         IC = M%CELL_INDEX(I-1,J,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE ! Permanently covered face
         IOR = -1
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            M%WALL(IW)%WALL_INDEX = 0
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I2
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (I==M%IBAR) CYCLE
         IC = M%CELL_INDEX(I+1,J,K)
         ! Permanently covered face
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE
         IOR = 1
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            M%WALL(IW)%WALL_INDEX = 0
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J1+1
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (J==1) CYCLE
         IC = M%CELL_INDEX(I,J-1,K)
         ! Permanently covered face
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE
         IOR = -2
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            M%WALL(IW)%WALL_INDEX = 0
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J2
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (J==M%JBAR) CYCLE
         IC = M%CELL_INDEX(I,J+1,K)
         ! Permanently covered face
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE
         IOR = 2
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            M%WALL(IW)%WALL_INDEX = 0
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K1+1
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (K==1) CYCLE
         IC = M%CELL_INDEX(I,J,K-1)
         ! Permanently covered face
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE
         IOR = -3
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            M%WALL(IW)%WALL_INDEX = 0
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K2
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (K==M%KBAR) CYCLE
         IC = M%CELL_INDEX(I,J,K+1)
         ! Permanently covered face
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE
         IOR = 3
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            M%WALL(IW)%WALL_INDEX = 0
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

ENDDO OBST_LOOP_2

! For MULT/OBST/SHAPE, recompute WC%ONE_D%AREA_ADJUST

OBST_SHAPE_IF: IF (OBST_SHAPE_AREA_ADJUST) THEN

   ! First, sum the face areas of the OBSTs with a given SURF
   SHAPE_LOOP_1: DO IW=M%N_EXTERNAL_WALL_CELLS+1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
      WC=>M%WALL(IW)
      IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE SHAPE_LOOP_1
      II=WC%ONE_D%II
      JJ=WC%ONE_D%JJ
      KK=WC%ONE_D%KK
      IOR=WC%ONE_D%IOR
      IC=M%CELL_INDEX(II,JJ,KK)
      I_OBST=M%OBST_INDEX_C(IC); IF (I_OBST==0)       CYCLE SHAPE_LOOP_1
      OB=>M%OBSTRUCTION(I_OBST); IF (OB%MULT_INDEX<0) CYCLE SHAPE_LOOP_1
      MR=>MULTIPLIER(OB%MULT_INDEX)
      SHAPE_SELECT_1: SELECT CASE(OB%SHAPE_TYPE)
         CASE(OBST_SPHERE_TYPE)
            MR%FDS_AREA(1) = MR%FDS_AREA(1) + WC%ONE_D%AREA
         CASE(OBST_CYLINDER_TYPE)
            ! OB%SHAPE_AREA follows the same pattern as SURF_IDS: top, sides, bottom
            SELECT CASE(IOR)
               CASE(3);         MR%FDS_AREA(1) = MR%FDS_AREA(1) + WC%ONE_D%AREA ! top
               CASE(-1,1,-2,2); MR%FDS_AREA(2) = MR%FDS_AREA(2) + WC%ONE_D%AREA ! side
               CASE(-3);        MR%FDS_AREA(3) = MR%FDS_AREA(3) + WC%ONE_D%AREA ! bottom
            END SELECT
         CASE(OBST_CONE_TYPE)
            SELECT CASE(IOR)
               CASE(-1,1,-2,2,3); MR%FDS_AREA(1) = MR%FDS_AREA(1) + WC%ONE_D%AREA
               CASE(-3);          MR%FDS_AREA(2) = MR%FDS_AREA(2) + WC%ONE_D%AREA
            END SELECT
         CASE(OBST_BOX_TYPE)
            ! Follows sextuplet ordering from SURF_ID6
            SELECT CASE(IOR)
               CASE(-1); MR%FDS_AREA(1) = MR%FDS_AREA(1) + WC%ONE_D%AREA
               CASE( 1); MR%FDS_AREA(2) = MR%FDS_AREA(2) + WC%ONE_D%AREA
               CASE(-2); MR%FDS_AREA(3) = MR%FDS_AREA(3) + WC%ONE_D%AREA
               CASE( 2); MR%FDS_AREA(4) = MR%FDS_AREA(4) + WC%ONE_D%AREA
               CASE(-3); MR%FDS_AREA(5) = MR%FDS_AREA(5) + WC%ONE_D%AREA
               CASE( 3); MR%FDS_AREA(6) = MR%FDS_AREA(6) + WC%ONE_D%AREA
            END SELECT
      END SELECT SHAPE_SELECT_1
   ENDDO SHAPE_LOOP_1

   ! Next, AREA_ADJUST the WALL_CELL
   SHAPE_LOOP_2: DO IW=M%N_EXTERNAL_WALL_CELLS+1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
      WC=>M%WALL(IW)
      IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE SHAPE_LOOP_2
      II=WC%ONE_D%II
      JJ=WC%ONE_D%JJ
      KK=WC%ONE_D%KK
      IOR=WC%ONE_D%IOR
      IC=M%CELL_INDEX(II,JJ,KK)
      I_OBST=M%OBST_INDEX_C(IC); IF (I_OBST==0)       CYCLE SHAPE_LOOP_2
      OB=>M%OBSTRUCTION(I_OBST); IF (OB%MULT_INDEX<0) CYCLE SHAPE_LOOP_2
      MR=>MULTIPLIER(OB%MULT_INDEX)
      SF=>SURFACE(WC%SURF_INDEX)
      IF (OB%SHAPE_TYPE>0) WC%ONE_D%AREA_ADJUST=1._EB
      SHAPE_SELECT_2: SELECT CASE(OB%SHAPE_TYPE)
         CASE(OBST_SPHERE_TYPE)
            WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1)
         CASE(OBST_CYLINDER_TYPE)
            SELECT CASE(IOR)
               CASE(3);         WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1) ! top
               CASE(-1,1,-2,2); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(2) ! side
               CASE(-3);        WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(3)/MR%FDS_AREA(3) ! bottom
            END SELECT
         CASE(OBST_CONE_TYPE)
            SELECT CASE(IOR)
               CASE(-1,1,-2,2,3); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1)
               CASE(-3);          WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(2)
            END SELECT
         CASE(OBST_BOX_TYPE)
            SELECT CASE(IOR)
               CASE(-1); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1)
               CASE( 1); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(2)
               CASE(-2); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(3)/MR%FDS_AREA(3)
               CASE( 2); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(4)
               CASE(-3); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(5)
               CASE( 3); WC%ONE_D%AREA_ADJUST = OB%SHAPE_AREA(3)/MR%FDS_AREA(6)
            END SELECT
      END SELECT SHAPE_SELECT_2
      WC%ONE_D%AREA_ADJUST = WC%ONE_D%AREA_ADJUST*SF%AREA_MULTIPLIER
   ENDDO SHAPE_LOOP_2

ENDIF OBST_SHAPE_IF

NOT_EVAC_IF_3: IF (.NOT.EVACUATION_ONLY(NM)) THEN
! Solid 3D heat and mass transfer

IF (ANY(OBSTRUCTION%HT3D)) SOLID_HT3D=.TRUE.
IF (ANY(OBSTRUCTION%MT3D)) SOLID_MT3D=.TRUE.

! Solid phase chemical heat source term

IF (SOLID_HT3D) THEN
   ALLOCATE(M%Q_DOT_PPP_S(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','Q_DOT_PPP_S',IZERO)
   M%Q_DOT_PPP_S = 0._EB
   ! Set internal heat source for 3d heat transfer
   OBST_LOOP_3: DO N=1,M%N_OBST
      OB=>M%OBSTRUCTION(N)
      DO K=OB%K1+1,OB%K2
         DO J=OB%J1+1,OB%J2
            DO I=OB%I1+1,OB%I2
               M%Q_DOT_PPP_S(I,J,K) = OB%INTERNAL_HEAT_SOURCE
            ENDDO
         ENDDO
      ENDDO
   ENDDO OBST_LOOP_3
ENDIF

! Assign TMP_INNER for HT3D (overrides INIT)

IF (SOLID_HT3D) THEN
   OBST_LOOP_4: DO N=1,M%N_OBST
      OB=>M%OBSTRUCTION(N)
      IF (OB%MATL_SURF_INDEX<=0) CYCLE OBST_LOOP_4
      SF=>SURFACE(OB%MATL_SURF_INDEX)
      DO K=OB%K1+1,OB%K2
         DO J=OB%J1+1,OB%J2
            DO I=OB%I1+1,OB%I2
               M%TMP(I,J,K) = SF%TMP_INNER(1)
            ENDDO
         ENDDO
      ENDDO
   ENDDO OBST_LOOP_4
ENDIF

! Solid 3D pyrolysis

IF (ANY(M%OBSTRUCTION%PYRO3D)) SOLID_PYRO3D=.TRUE.

! Initialize material densities for 3D pyrolysis

IF (SOLID_PYRO3D) THEN
   OBST_LOOP_5: DO N=1,M%N_OBST
      OB=>M%OBSTRUCTION(N)
      IF (OB%MATL_SURF_INDEX<=0) CYCLE OBST_LOOP_5
      SF=>SURFACE(OB%MATL_SURF_INDEX)
      MATL_LOOP: DO NN = 1,SF%N_MATL
         DO K=OB%K1+1,OB%K2
            DO J=OB%J1+1,OB%J2
               DO I=OB%I1+1,OB%I2
                  OB%RHO(I,J,K,NN) = SF%RHO_0(1,NN)
                  IF (MATERIAL(NN)%PYROLYSIS_MODEL==PYROLYSIS_LIQUID) OB%PYRO3D_LIQUID=.TRUE.
               ENDDO
            ENDDO
         ENDDO
      ENDDO MATL_LOOP
   ENDDO OBST_LOOP_5
ENDIF

! Initialize gas densities for 3D pyrolysis

IF (SOLID_MT3D) THEN
   ! store mass production rate for tracked species
   ALLOCATE(M%M_DOT_G_PPP_S(0:M%IBP1,0:M%JBP1,0:M%KBP1,1:N_TRACKED_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','M_DOT_G_PPP_S',IZERO)
   M%M_DOT_G_PPP_S = 0._EB
   ! initialize gas phase mass densities
   ALLOCATE(M%RHO_ZZ_G_S(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','RHO_ZZ_G_S',IZERO)
   M%RHO_ZZ_G_S = 0._EB
   ! !! debug !!
   ! ALLOCATE(M%SCALAR_SAVE1(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES),STAT=IZERO)
   ! CALL ChkMemErr('INIT','SCALAR_SAVE1',IZERO)
   ! M%SCALAR_SAVE1 = 0._EB
ENDIF

! Store old velocity components for special output quantities

IF (STORE_OLD_VELOCITY) THEN
   ALLOCATE(M%U_OLD(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO); CALL ChkMemErr('INIT','U_OLD',IZERO)
   ALLOCATE(M%V_OLD(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO); CALL ChkMemErr('INIT','V_OLD',IZERO)
   ALLOCATE(M%W_OLD(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO); CALL ChkMemErr('INIT','W_OLD',IZERO)
ENDIF

! Wind

IF (ANY(MEAN_FORCING)) THEN
   DO IW=1,M%N_EXTERNAL_WALL_CELLS
      WC => M%WALL(IW)
      IOR = WC%ONE_D%IOR
      IF (SPONGE_CELLS>0 .AND. WC%BOUNDARY_TYPE==OPEN_BOUNDARY .AND. MEAN_FORCING(ABS(IOR))) M%APPLY_SPONGE_LAYER(IOR) = .TRUE.
   ENDDO
ENDIF

! Reset ghost cell values of cell centered velocity for use in computing viscosity (must be done after INIT_WALL_CELL)

DO IW=1,M%N_EXTERNAL_WALL_CELLS
   M%U_GHOST(IW) = U0*EVALUATE_RAMP(T_BEGIN,DUMMY,I_RAMP_U0_T)*EVALUATE_RAMP(M%ZC(M%WALL(IW)%ONE_D%KKG),DUMMY,I_RAMP_U0_Z)
   M%V_GHOST(IW) = V0*EVALUATE_RAMP(T_BEGIN,DUMMY,I_RAMP_V0_T)*EVALUATE_RAMP(M%ZC(M%WALL(IW)%ONE_D%KKG),DUMMY,I_RAMP_V0_Z)
   M%W_GHOST(IW) = W0*EVALUATE_RAMP(T_BEGIN,DUMMY,I_RAMP_W0_T)*EVALUATE_RAMP(M%ZC(M%WALL(IW)%ONE_D%KKG),DUMMY,I_RAMP_W0_Z)
ENDDO

! Find exterior wall cells of HT3D solids

IF (SOLID_PYRO3D .AND. .NOT.SOLID_MT3D) CALL ASSIGN_HT3D_WALL_INDICES(NM)

ENDIF NOT_EVAC_IF_3

END SUBROUTINE INITIALIZE_MESH_VARIABLES_1


SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM)

USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_HEAT,LES_FILTER_WIDTH_FUNCTION
USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
USE CONTROL_VARIABLES
INTEGER :: N,I,J,K,II,JJ,KK,IPTS,JPTS,KPTS,N_EDGES_DIM,IW,IC,ICG,IOR,IERR,IPZ,NOM,ITER,IZERO,IIG,JJG,KKG,ICF,NSLICE
INTEGER, INTENT(IN) :: NM
REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),VC,RTRM,CP,XXC,YYC,ZZC
INTEGER, POINTER :: IBP1, JBP1, KBP1,IBAR, JBAR, KBAR, N_EDGES
INTEGER, ALLOCATABLE, DIMENSION(:) :: IW_EXPORT
REAL(EB),POINTER :: XS,XF,YS,YF,ZS,ZF
TYPE (MESH_TYPE), POINTER :: M
TYPE (WALL_TYPE), POINTER :: WC
TYPE (CFACE_TYPE), POINTER :: CFA
TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE (MESH_TYPE), POINTER :: OM
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (SURFACE_TYPE), POINTER :: SF
LOGICAL :: SOLID_CELL

IERR = 0
M => MESHES(NM)
IBP1 =>M%IBP1
JBP1 =>M%JBP1
KBP1 =>M%KBP1
IBAR =>M%IBAR
JBAR =>M%JBAR
KBAR =>M%KBAR
N_EDGES=>M%N_EDGES
XS=>M%XS
YS=>M%YS
ZS=>M%ZS
XF=>M%XF
YF=>M%YF
ZF=>M%ZF

! Surface work arrays

ALLOCATE(M%WALL_WORK1(M%N_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK1',IZERO)
ALLOCATE(M%WALL_WORK2(M%N_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK2',IZERO)

! Background pressure variables

ALLOCATE(  M%PBAR(0:KBP1,0:N_ZONE),STAT=IZERO)   ; CALL ChkMemErr('INIT','PBAR',IZERO)
ALLOCATE(  M%PBAR_S(0:KBP1,0:N_ZONE),STAT=IZERO) ; CALL ChkMemErr('INIT','PBAR_S',IZERO)
ALLOCATE(  M%R_PBAR(0:KBP1,0:N_ZONE),STAT=IZERO) ; CALL ChkMemErr('INIT','R_PBAR',IZERO)
ALLOCATE(  M%D_PBAR_DT(0:N_ZONE),STAT=IZERO)     ; CALL ChkMemErr('INIT','D_PBAR_DT',IZERO)    ; M%D_PBAR_DT = 0._EB
ALLOCATE(  M%D_PBAR_DT_S(0:N_ZONE),STAT=IZERO)   ; CALL ChkMemErr('INIT','D_PBAR_DT_S',IZERO)  ; M%D_PBAR_DT_S = 0._EB
ALLOCATE(  M%U_LEAK(0:N_ZONE),STAT=IZERO)        ; CALL ChkMemErr('INIT','U_LEAK',IZERO)       ; M%U_LEAK = 0._EB

DO K=0,M%KBP1
   M%PBAR(K,:)   = M%P_0(K)
   M%PBAR_S(K,:) = M%P_0(K)
ENDDO

! Initialize PSUM for zone cases

IF (N_ZONE > 0) THEN
   N_ZONE_LOOP: DO IPZ = 1,N_ZONE
      PSUM(IPZ,NM) = 0._EB
      IF (EVACUATION_ONLY(NM)) EXIT N_ZONE_LOOP
      DO K=1,M%KBAR
         DO J=1,M%JBAR
            DO I=1,M%IBAR
               IF (M%PRESSURE_ZONE(I,J,K) /= IPZ) CYCLE
               IF (M%SOLID(M%CELL_INDEX(I,J,K)))    CYCLE
               VC   = M%DX(I)*M%RC(I)*M%DY(J)*M%DZ(K)
               ZZ_GET(1:N_TRACKED_SPECIES) = M%ZZ(I,J,K,1:N_TRACKED_SPECIES)
               CALL GET_SPECIFIC_HEAT(ZZ_GET,CP,M%TMP(I,J,K))
               RTRM = M%RSUM(I,J,K)/(CP*M%PBAR(K,IPZ))
               PSUM(IPZ,NM) = PSUM(IPZ,NM) + VC*(1._EB/M%PBAR(K,IPZ)-RTRM)
            ENDDO
         ENDDO
      ENDDO
   ENDDO N_ZONE_LOOP
ENDIF

! Loop through WALL and CFACE cells and assign PRESSURE_ZONE. Also, check for
! inappropriate boundaries, such as thin obstructions that burn or blow.

WALL_LOOP_0: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS

   WC => M%WALL(IW)
   SF => SURFACE(WC%SURF_INDEX)

   IF (.NOT.EVACUATION_ONLY(NM)) WC%ONE_D%PRESSURE_ZONE = M%PRESSURE_ZONE(WC%ONE_D%IIG,WC%ONE_D%JJG,WC%ONE_D%KKG)

   IF (IW<=M%N_EXTERNAL_WALL_CELLS) THEN
      IF (M%EXTERNAL_WALL(IW)%NOM>0) THEN
         OM => MESHES(M%EXTERNAL_WALL(IW)%NOM)
         IC = OM%CELL_INDEX(M%EXTERNAL_WALL(IW)%IIO_MIN,M%EXTERNAL_WALL(IW)%JJO_MIN,M%EXTERNAL_WALL(IW)%KKO_MIN)
         SOLID_CELL = OM%SOLID(IC)
      ELSE
         IC =  M%CELL_INDEX(WC%ONE_D%II,WC%ONE_D%JJ,WC%ONE_D%KK)
         SOLID_CELL = M%SOLID(IC)
      ENDIF
   ELSE
      IC =  M%CELL_INDEX(WC%ONE_D%II,WC%ONE_D%JJ,WC%ONE_D%KK)
      SOLID_CELL = M%SOLID(IC)
   ENDIF

   IF (.NOT.SOLID_CELL) THEN
      IF ( (ABS(WC%ONE_D%U_NORMAL_0)>TWO_EPSILON_EB .OR. ANY(SF%LEAK_PATH>=0)) .AND. WC%OBST_INDEX>0 ) THEN
         WRITE(LU_ERR,'(A,A,A,I0)') 'ERROR: SURF ',TRIM(SF%ID),' cannot be applied to a thin obstruction, OBST #',&
                                    M%OBSTRUCTION(WC%OBST_INDEX)%ORDINAL
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF
      IF (WC%VENT_INDEX>0 .AND. WC%OBST_INDEX>0) THEN
         VT => VENTS(WC%VENT_INDEX)
         IF (VT%BOUNDARY_TYPE==HVAC_BOUNDARY) THEN
            WRITE(LU_ERR,'(A,A,A,I0)') 'ERROR: VENT ',TRIM(VT%ID),' cannot be applied to a thin obstruction, OBST #',&
                                    M%OBSTRUCTION(WC%OBST_INDEX)%ORDINAL
            STOP_STATUS = SETUP_STOP
            RETURN
         ENDIF
      ENDIF
   ENDIF

ENDDO WALL_LOOP_0

CFACE_LOOP_0: DO ICF=1,M%N_CFACE_CELLS
   CFA => M%CFACE(ICF)
   CFA%ONE_D%PRESSURE_ZONE = M%PRESSURE_ZONE(CFA%ONE_D%IIG,CFA%ONE_D%JJG,CFA%ONE_D%KKG)
ENDDO CFACE_LOOP_0

! If there is complex terrain using GEOM and Above Ground Level (AGL) slices,
! determine K index of gas phase quantities.

IF (CC_IBM) THEN
   DO NSLICE = 1, M%N_TERRAIN_SLCF
      DO ICF=1,M%N_CUTFACE_MESH
         IF (CUT_FACE(ICF)%STATUS/=2 .OR. CUT_FACE(ICF)%NFACE<1) CYCLE
         IW  = MAXLOC(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE),DIM=1)
         CFA => CFACE( CUT_FACE(ICF)%CFACE_INDEX(IW) )
         IF (CFA%NVEC(KAXIS)>-TWO_EPSILON_EB .AND.  CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
            IF (CFA%ONE_D%KKG > M%K_AGL_SLICE(CFA%ONE_D%IIG,CFA%ONE_D%JJG,NSLICE)) THEN
               M%K_AGL_SLICE(CFA%ONE_D%IIG,CFA%ONE_D%JJG,NSLICE) = M%K_AGL_SLICE(CFA%ONE_D%IIG,CFA%ONE_D%JJG,NSLICE) + CFA%ONE_D%KKG
            ENDIF
         ENDIF
      ENDDO
   ENDDO
ENDIF

! Loop through all internal and external wall cells and look for thermally thick
! solids with EXPOSED back wall cells. If the exposed back wall cell is in
! another mesh, store the cell info into arrays that are to be MPI exchanged.

! DEFINITION MESHES(NM)%OMESH(NOM)%N_EXPOSED_WALL_CELLS
! Number of wall cells in Mesh NM whose exposed back faces are in Mesh NOM.

! DEFINITION MESHES(NM)%OMESH(NOM)%EXPOSED_WALL_CELL_BACK_INDICES(1:N_EXPOSED_WALL_CELLS)
! Wall indices of the back faces of the exposed wall cells.

NON_EVAC_IF: IF (.NOT.EVACUATION_ONLY(NM)) THEN

DO NOM=1,NMESHES
   M%OMESH(NOM)%N_EXPOSED_WALL_CELLS = 0
ENDDO

WALL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS

   WC => M%WALL(IW)

   II  = WC%ONE_D%II
   JJ  = WC%ONE_D%JJ
   KK  = WC%ONE_D%KK
   IIG = WC%ONE_D%IIG
   JJG = WC%ONE_D%JJG
   KKG = WC%ONE_D%KKG
   ICG = M%CELL_INDEX(IIG,JJG,KKG)
   IOR = WC%ONE_D%IOR

   ! Special case where an HT3D solid abuts a mesh boundary and needs temperature and heat flux from the other mesh.

   IF (M%OBST_INDEX_C(ICG)/=0 .AND. M%OBSTRUCTION(M%OBST_INDEX_C(ICG))%HT3D .AND. IW<=M%N_EXTERNAL_WALL_CELLS) THEN
      EWC => M%EXTERNAL_WALL(IW)
      NOM = EWC%NOM
      IF (NOM/=NM .AND. NOM>0) THEN
         M%OMESH(NOM)%N_EXPOSED_WALL_CELLS = M%OMESH(NOM)%N_EXPOSED_WALL_CELLS + 1
         IC = MESHES(NOM)%CELL_INDEX(EWC%IIO_MIN,EWC%JJO_MIN,EWC%KKO_MIN)
         WC%BACK_INDEX = MESHES(NOM)%WALL_INDEX(IC,IOR)
         WC%BACK_MESH  = NOM
         CYCLE WALL_LOOP
      ENDIF
   ENDIF

   ! Search for 0 or 1 cell thick HT1D solids that abut mesh boundary and have EXPOSED back boundary condition.

   IF (WC%OBST_INDEX==0) CYCLE WALL_LOOP
   IF (M%SOLID(ICG)) CYCLE WALL_LOOP
   SF => SURFACE(WC%SURF_INDEX)

   IF_THERM_THICK_EXPOSED: IF (SF%THERMAL_BC_INDEX==THERMALLY_THICK .AND. SF%BACKING==EXPOSED) THEN ! search for the back side

      NOM = NM
      OM => MESHES(NOM)
      II = WC%ONE_D%II
      JJ = WC%ONE_D%JJ
      KK = WC%ONE_D%KK

      ZERO_OR_ONE_CELL_THICK: DO ITER=1,2  ! Look for the back wall cell face if the obstruction is either zero or one cell thick.

         IF (II==0 .OR. II==OM%IBP1 .OR. JJ==0 .OR. JJ==OM%JBP1 .OR. KK==0 .OR. KK==OM%KBP1) THEN
            XXC=OM%XC(II) ; YYC=OM%YC(JJ) ; ZZC=OM%ZC(KK)
            IF (II==0)       XXC = OM%X(II)   - MESH_SEPARATION_DISTANCE
            IF (II==OM%IBP1) XXC = OM%X(II-1) + MESH_SEPARATION_DISTANCE
            IF (JJ==0)       YYC = OM%Y(JJ)   - MESH_SEPARATION_DISTANCE
            IF (JJ==OM%JBP1) YYC = OM%Y(JJ-1) + MESH_SEPARATION_DISTANCE
            IF (KK==0)       ZZC = OM%Z(KK)   - MESH_SEPARATION_DISTANCE
            IF (KK==OM%KBP1) ZZC = OM%Z(KK-1) + MESH_SEPARATION_DISTANCE
            CALL SEARCH_OTHER_MESHES(XXC,YYC,ZZC,NOM,II,JJ,KK)
            IF (NOM==0) CYCLE WALL_LOOP
            OM => MESHES(NOM)
         ENDIF
         IC = OM%CELL_INDEX(II,JJ,KK)
         IF (.NOT.OM%SOLID(IC) .AND. OM%WALL_INDEX(IC,IOR)>0) THEN ! the back wall face is found
            WC%BACK_INDEX = OM%WALL_INDEX(IC,IOR)
            WC%BACK_MESH  = NOM
            EXIT ZERO_OR_ONE_CELL_THICK
         ELSE ! see if the obstruction is one cell thick. If not, get out.
            IF (ITER==2) CYCLE WALL_LOOP
            SELECT CASE(IOR)
               CASE(-1) ; II=II+1
               CASE( 1) ; II=II-1
               CASE(-2) ; JJ=JJ+1
               CASE( 2) ; JJ=JJ-1
               CASE(-3) ; KK=KK+1
               CASE( 3) ; KK=KK-1
            END SELECT
         ENDIF

      ENDDO ZERO_OR_ONE_CELL_THICK

      IF (NOM/=NM) THEN ! count the number of wall cells that need to be exchanged
         M%OMESH(NOM)%N_EXPOSED_WALL_CELLS = M%OMESH(NOM)%N_EXPOSED_WALL_CELLS + 1
      ENDIF

   ENDIF IF_THERM_THICK_EXPOSED

ENDDO WALL_LOOP

DO NOM=1,NMESHES
   IF (M%OMESH(NOM)%N_EXPOSED_WALL_CELLS==0) CYCLE
   ALLOCATE(M%OMESH(NOM)%EXPOSED_WALL_CELL_BACK_INDICES(M%OMESH(NOM)%N_EXPOSED_WALL_CELLS))
ENDDO

ALLOCATE(IW_EXPORT(NMESHES)) ; IW_EXPORT = 0  ! This is just a counter for each neighboring mesh.

DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
   WC => M%WALL(IW)
   IF (WC%BACK_MESH/=NM) THEN ! Save the back wall cell index, and then reset the index to form a short list of exposed back cells.
      IW_EXPORT(WC%BACK_MESH) = IW_EXPORT(WC%BACK_MESH) + 1
      M%OMESH(WC%BACK_MESH)%EXPOSED_WALL_CELL_BACK_INDICES(IW_EXPORT(WC%BACK_MESH)) = WC%BACK_INDEX
      WC%BACK_INDEX = IW_EXPORT(WC%BACK_MESH)
   ENDIF
ENDDO

DEALLOCATE(IW_EXPORT)

ENDIF NON_EVAC_IF

! Set clocks and counters related to frequency of solid phase conduction updates

M%BC_CLOCK     = T_BEGIN
M%WALL_COUNTER = 0
M%BC_CLOCK_HT3D     = T_BEGIN
M%WALL_COUNTER_HT3D = 0

! Allocate arrays for storing velocity boundary condition info

N_EDGES_DIM = 4*(IBP1*JBP1+IBP1*KBP1+JBP1*KBP1)
IF (EVACUATION_ONLY(NM)) N_EDGES_DIM = 4*(IBP1*KBP1+JBP1*KBP1)
DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
   IPTS = OB%I2-OB%I1
   JPTS = OB%J2-OB%J1
   KPTS = OB%K2-OB%K1
   IF (EVACUATION_ONLY(NM)) THEN
      N_EDGES_DIM = N_EDGES_DIM + 4*(IPTS*KPTS+JPTS*KPTS)
   ELSE
      N_EDGES_DIM = N_EDGES_DIM + 4*(IPTS*JPTS+IPTS*KPTS+JPTS*KPTS)
   ENDIF
ENDDO

ALLOCATE(M%IJKE(16,N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','IJKE',IZERO)
M%IJKE  = 0
ALLOCATE(M%OME_E(-2:2,0:N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','OME_E',IZERO)
M%OME_E = 0._EB
ALLOCATE(M%TAU_E(-2:2,0:N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','TAU_E',IZERO)
M%TAU_E = 0._EB
ALLOCATE(M%EDGE_INTERPOLATION_FACTOR(N_EDGES_DIM,2),STAT=IZERO)
CALL ChkMemErr('INIT','EDGE_INTERPOLATION_FACTOR',IZERO)
M%EDGE_INTERPOLATION_FACTOR = 1._EB

! Allocate array to hold character strings for Smokeview file

M%N_STRINGS     =   0
M%N_STRINGS_MAX = 100
ALLOCATE(M%STRING(M%N_STRINGS_MAX),STAT=IZERO)
CALL ChkMemErr('INIT','STRING',IZERO)

! Set up arrays to hold velocity boundary condition info

CALL INITIALIZE_EDGES

! Initialize Pressure solver

IF (.NOT.FREEZE_VELOCITY) CALL INITIALIZE_POISSON_SOLVER
!CALL INITIALIZE_LAPLACE_SOLVER
IF (IERR/=0) RETURN

! Initialize Mesh Exchange

CALL INITIALIZE_INTERPOLATION


CONTAINS


SUBROUTINE INITIALIZE_EDGES

! Set up edge arrays for velocity boundary conditions

INTEGER I,J,K,N

CALL POINT_TO_MESH(NM)

N_EDGES = 0

! Arguments for DEFINE_EDGE(I,J,K,IOR,IEC,NM,OBST_INDEX)

DO K=0,KBAR
   DO J=0,JBAR
      IF (J>0) CALL DEFINE_EDGE(   0,J,K, 1,2,NM,0,IERR)
      IF (J>0) CALL DEFINE_EDGE(IBAR,J,K,-1,2,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(   0,J,K, 1,3,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(IBAR,J,K,-1,3,NM,0,IERR)
   ENDDO
ENDDO
DO K=0,KBAR
   DO I=0,IBAR
      IF (I>0) CALL DEFINE_EDGE(I,   0,K, 2,1,NM,0,IERR)
      IF (I>0) CALL DEFINE_EDGE(I,JBAR,K,-2,1,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(I,   0,K, 2,3,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(I,JBAR,K,-2,3,NM,0,IERR)
   ENDDO
ENDDO
DO J=0,JBAR
   DO I=0,IBAR
      IF (I>0) CALL DEFINE_EDGE(I,J,   0, 3,1,NM,0,IERR)
      IF (I>0) CALL DEFINE_EDGE(I,J,KBAR,-3,1,NM,0,IERR)
      IF (J>0) CALL DEFINE_EDGE(I,J,   0, 3,2,NM,0,IERR)
      IF (J>0) CALL DEFINE_EDGE(I,J,KBAR,-3,2,NM,0,IERR)
   ENDDO
ENDDO

IF (IERR/=0) RETURN

OBST_LOOP_3: DO N=1,N_OBST
   OB => OBSTRUCTION(N)
   DO K=OB%K1,OB%K2
      DO J=OB%J1,OB%J2
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I1,J,K,-1,2,NM,N,IERR)
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I2,J,K, 1,2,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I1,J,K,-1,3,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I2,J,K, 1,3,NM,N,IERR)
      ENDDO
   ENDDO
   DO K=OB%K1,OB%K2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J1,K,-2,1,NM,N,IERR)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J2,K, 2,1,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J1,K,-2,3,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J2,K, 2,3,NM,N,IERR)
      ENDDO
   ENDDO
   DO J=OB%J1,OB%J2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K1,-3,1,NM,N,IERR)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K2, 3,1,NM,N,IERR)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K1,-3,2,NM,N,IERR)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K2, 3,2,NM,N,IERR)
      ENDDO
   ENDDO
ENDDO OBST_LOOP_3

IF (N_EDGES>N_EDGES_DIM .AND. EVACUATION_ONLY(NM)) THEN
   WRITE(LU_ERR,'(A,I2,A,2I8)') 'ERROR: Edges memory; Mesh: ',NM,', n_edges, n_edges_dim ',N_EDGES, N_EDGES_DIM
   STOP_STATUS = SETUP_STOP
   IERR = 1
ENDIF

END SUBROUTINE INITIALIZE_EDGES


SUBROUTINE INITIALIZE_POISSON_SOLVER

USE POIS, ONLY: H3CZIS,H2CZIS,H3CSIS,H2CYIS
REAL(EB) :: XLM,XMU
INTEGER  :: N,IZERO
INTEGER, POINTER :: ITRN,JTRN,KTRN,LBC,MBC,NBC
INTEGER, POINTER, DIMENSION(:) :: NOC
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (WALL_TYPE), POINTER :: WC

! Allocate major arrays

ITRN =>M%ITRN
JTRN =>M%JTRN
KTRN =>M%KTRN
LBC =>M%LBC
MBC =>M%MBC
NBC =>M%NBC
LBC2 =>M%LBC2
MBC2 =>M%MBC2
NBC2 =>M%NBC2
NOC=>TRANS(NM)%NOC
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=0
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=1
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=2
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=3
IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=4
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=5
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) M%IPS=6
IF (EVACUATION_ONLY(NM)                      ) M%IPS=7
IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) THEN
   WRITE(LU_ERR,'(A,I3)') 'ERROR: Stretch at most 2 coordinate directions in MESH ',NM
   STOP_STATUS = SETUP_STOP
   IERR = 1
   RETURN
ENDIF

IF (M%IPS<=1 .OR. M%IPS==4) THEN
   ITRN = IBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = KBP1

   ! pressure periodic boundary conditions
   IF (FISHPAK_BC(1)==FISHPAK_BC_PERIODIC) ITRN=IBAR
   IF (FISHPAK_BC(2)==FISHPAK_BC_PERIODIC) JTRN=JBAR
   IF (FISHPAK_BC(3)==FISHPAK_BC_PERIODIC) KTRN=KBAR
ENDIF

IF (M%IPS==2) THEN
   ITRN = JBP1
   JTRN = IBP1
   KTRN = KBP1
   ALLOCATE(M%BZST(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF

IF (M%IPS==3 .OR. M%IPS==6) THEN
   ITRN = KBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = IBP1
   ALLOCATE(M%BXST(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
   ALLOCATE(M%BYST(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYST',IZERO)
   ALLOCATE(M%BYFT(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYFT',IZERO)
   ALLOCATE(M%BZST(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF

IF (M%IPS==5) THEN
   ITRN = IBP1
   JTRN = KBP1
   KTRN = JBP1
   ALLOCATE(M%BXST(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
ENDIF

IF (M%IPS==7) THEN
   ITRN = IBP1
   JTRN = JBP1
   KTRN = 1
ENDIF

IF (M%IPS<=3 .OR. M%IPS==7) THEN
   M%LSAVE = (ITRN+1)*JTRN*KTRN+7*ITRN+5*JTRN+6*KTRN+56
   M%LWORK = (ITRN+1)*JTRN*KTRN
ELSE
   N_LOOP: DO N=1,50
      IF ((JTRN+1)<=2**N) EXIT N_LOOP
   ENDDO N_LOOP
   M%LSAVE = KTRN*(6*N*(2**N)+2*N+19)+8*ITRN+7*JTRN+38
   M%LWORK = JTRN*(ITRN*(KTRN+1)+1)
ENDIF

ALLOCATE(M%SAVE1(-3:M%LSAVE),STAT=IZERO)
CALL ChkMemErr('INIT','SAVE1',IZERO)
IF (EXTERNAL_BOUNDARY_CORRECTION) THEN
   ALLOCATE(M%SAVE2(-3:M%LSAVE),STAT=IZERO)
   CALL ChkMemErr('INIT','SAVE2',IZERO)
ENDIF
ALLOCATE(M%WORK(M%LWORK),STAT=IZERO)
CALL ChkMemErr('INIT','WORK',IZERO)
ALLOCATE(M%PRHS(ITRN,JTRN,KTRN),STAT=IZERO)
CALL ChkMemErr('INIT','PRHS',IZERO)
IF (KBAR>1) THEN
   IF (JBAR>1 ) ALLOCATE(M%BXS(JBP1,KBP1),STAT=IZERO)
   IF (JBAR==1) ALLOCATE(M%BXS(1,KBP1)   ,STAT=IZERO)
ELSE
   ALLOCATE(M%BXS(JBP1,1)   ,STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BXS',IZERO)
IF (KBAR>1) THEN
   IF (JBAR>1 ) ALLOCATE(M%BXF(JBP1,KBP1),STAT=IZERO)
   IF (JBAR==1) ALLOCATE(M%BXF(1,KBP1)   ,STAT=IZERO)
ELSE
   ALLOCATE(M%BXF(JBP1,1)   ,STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BXF',IZERO)
IF (KBAR>1) THEN
   ALLOCATE(M%BYS(IBP1,KBP1),STAT=IZERO)
ELSE
   ALLOCATE(M%BYS(IBP1,1),STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BYS',IZERO)
IF (KBAR>1) THEN
   ALLOCATE(M%BYF(IBP1,KBP1),STAT=IZERO)
ELSE
   ALLOCATE(M%BYF(IBP1,1),STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BYF',IZERO)
IF (JBAR>1) ALLOCATE(M%BZS(IBP1,JBP1),STAT=IZERO)
IF (JBAR==1) ALLOCATE(M%BZS(IBP1,1)   ,STAT=IZERO)
CALL ChkMemErr('INIT','BZS',IZERO)
IF (JBAR>1) ALLOCATE(M%BZF(IBP1,JBP1),STAT=IZERO)
IF (JBAR==1) ALLOCATE(M%BZF(IBP1,1)   ,STAT=IZERO)
CALL ChkMemErr('INIT','BZF',IZERO)

M%POIS_PTB = 0._EB
M%SAVE1 = 0._EB
M%WORK = 0._EB
M%PRHS = 0._EB
M%BXS  = 0._EB
M%BXF  = 0._EB
M%BYS  = 0._EB
M%BYF  = 0._EB
M%BZS  = 0._EB
M%BZF  = 0._EB

IF (EXTERNAL_BOUNDARY_CORRECTION) THEN
   M%LAPLACE_PTB = 0._EB
   M%SAVE2 = 0._EB
ENDIF

! Initialize pressure solver

XLM = 0._EB         ! No Helmholtz equation
XMU = 0._EB         ! No Helmholtz equation

! BC index for Fishpak solver

! From Roland Sweet's notes:
!
! Here we use LBC as an example, this is the BC index for the X direction. MBC and NMC are
! analogous for the Y and Z directions.
!
! LBC = 0, solution is periodic in X.
!
! LBC = 1, solution is specified at XS (Dirichlet) and XF (Dirichlet).
!
! LBC = 2, solution is specified at XS (Dirichlet) and derivative of solution is specified at XF (Neumann).
!
! LBC = 3, derivative of solution is specified at XS (Neumann) and XF (Neumann).
!
! LBC = 4, derivative of solution is specified at XS (Neumann) and solution is specified at XF (Dirichlet).
!
! LBC = 5, the solution is unspecified at r = RS = 0 and the solution is specified at r = RF.
!
! LBC = 6, if the solution is unspecified at r = RS = 0 and the derivative of the solution with respect to r is specified
!          at r = RF.

LBC = FISHPAK_BC_NEUMANN_NEUMANN
MBC = FISHPAK_BC_NEUMANN_NEUMANN
NBC = FISHPAK_BC_NEUMANN_NEUMANN
IF (EXTERNAL_BOUNDARY_CORRECTION) THEN
   LBC2 = FISHPAK_BC_NEUMANN_NEUMANN
   MBC2 = FISHPAK_BC_NEUMANN_NEUMANN
   NBC2 = FISHPAK_BC_NEUMANN_NEUMANN
ENDIF

! Look for OPEN vents -- this will change the entire face to DIRICHLET BCs

VENT_LOOP: DO N=1,M%N_VENT
   VT => M%VENTS(N)
   IF (VT%BOUNDARY_TYPE /= OPEN_BOUNDARY .OR. VT%GHOST_CELLS_ONLY) CYCLE VENT_LOOP
   IF (VT%I1==0 .AND. VT%I2==0) THEN
      IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_DIRICHLET_NEUMANN
      IF (LBC==FISHPAK_BC_NEUMANN_DIRICHLET) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%I1==M%IBAR .AND. VT%I2==M%IBAR) THEN
      IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_NEUMANN_DIRICHLET
      IF (LBC==FISHPAK_BC_DIRICHLET_NEUMANN) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%J1==0 .AND. VT%J2==0) THEN
      IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_DIRICHLET_NEUMANN
      IF (MBC==FISHPAK_BC_NEUMANN_DIRICHLET) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%J1==M%JBAR .AND. VT%J2==M%JBAR) THEN
      IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_NEUMANN_DIRICHLET
      IF (MBC==FISHPAK_BC_DIRICHLET_NEUMANN) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%K1==0 .AND. VT%K2==0) THEN
      IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_DIRICHLET_NEUMANN
      IF (NBC==FISHPAK_BC_NEUMANN_DIRICHLET) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%K1==M%KBAR .AND. VT%K2==M%KBAR) THEN
      IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_NEUMANN_DIRICHLET
      IF (NBC==FISHPAK_BC_DIRICHLET_NEUMANN) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
ENDDO VENT_LOOP

! All interpolated boundaries are Dirichlet

DO IW=1,M%N_EXTERNAL_WALL_CELLS
   IF (M%EXTERNAL_WALL(IW)%NOM==0) CYCLE
   SELECT CASE(M%WALL(IW)%ONE_D%IOR)
      CASE( 1)
         IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_DIRICHLET_NEUMANN
         IF (LBC==FISHPAK_BC_NEUMANN_DIRICHLET) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE(-1)
         IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_NEUMANN_DIRICHLET
         IF (LBC==FISHPAK_BC_DIRICHLET_NEUMANN) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE( 2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_DIRICHLET_NEUMANN
         IF (MBC==FISHPAK_BC_NEUMANN_DIRICHLET) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE(-2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_NEUMANN_DIRICHLET
         IF (MBC==FISHPAK_BC_DIRICHLET_NEUMANN) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE( 3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_DIRICHLET_NEUMANN
         IF (NBC==FISHPAK_BC_NEUMANN_DIRICHLET) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE(-3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_NEUMANN_DIRICHLET
         IF (NBC==FISHPAK_BC_DIRICHLET_NEUMANN) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   END SELECT
ENDDO

! Test setting interpolated boundaries to Dirichlet for Laplace solve

IF (EXTERNAL_BOUNDARY_CORRECTION) THEN
   DO IW=1,M%N_EXTERNAL_WALL_CELLS
      IF (M%EXTERNAL_WALL(IW)%NOM==0) CYCLE
      SELECT CASE(M%WALL(IW)%ONE_D%IOR)
         CASE( 1)
            IF (LBC2==FISHPAK_BC_NEUMANN_NEUMANN)   LBC2 = FISHPAK_BC_DIRICHLET_NEUMANN
            IF (LBC2==FISHPAK_BC_NEUMANN_DIRICHLET) LBC2 = FISHPAK_BC_DIRICHLET_DIRICHLET
         CASE(-1)
            IF (LBC2==FISHPAK_BC_NEUMANN_NEUMANN)   LBC2 = FISHPAK_BC_NEUMANN_DIRICHLET
            IF (LBC2==FISHPAK_BC_DIRICHLET_NEUMANN) LBC2 = FISHPAK_BC_DIRICHLET_DIRICHLET
         CASE( 2)
            IF (MBC2==FISHPAK_BC_NEUMANN_NEUMANN)   MBC2 = FISHPAK_BC_DIRICHLET_NEUMANN
            IF (MBC2==FISHPAK_BC_NEUMANN_DIRICHLET) MBC2 = FISHPAK_BC_DIRICHLET_DIRICHLET
         CASE(-2)
            IF (MBC2==FISHPAK_BC_NEUMANN_NEUMANN)   MBC2 = FISHPAK_BC_NEUMANN_DIRICHLET
            IF (MBC2==FISHPAK_BC_DIRICHLET_NEUMANN) MBC2 = FISHPAK_BC_DIRICHLET_DIRICHLET
         CASE( 3)
            IF (NBC2==FISHPAK_BC_NEUMANN_NEUMANN)   NBC2 = FISHPAK_BC_DIRICHLET_NEUMANN
            IF (NBC2==FISHPAK_BC_NEUMANN_DIRICHLET) NBC2 = FISHPAK_BC_DIRICHLET_DIRICHLET
         CASE(-3)
            IF (NBC2==FISHPAK_BC_NEUMANN_NEUMANN)   NBC2 = FISHPAK_BC_NEUMANN_DIRICHLET
            IF (NBC2==FISHPAK_BC_DIRICHLET_NEUMANN) NBC2 = FISHPAK_BC_DIRICHLET_DIRICHLET
      END SELECT
   ENDDO
ENDIF

! Periodic pressure boundary conditions for CrayFishpak

IF (FISHPAK_BC(1)==FISHPAK_BC_PERIODIC) LBC=FISHPAK_BC_PERIODIC
IF (FISHPAK_BC(2)==FISHPAK_BC_PERIODIC) MBC=FISHPAK_BC_PERIODIC
IF (FISHPAK_BC(3)==FISHPAK_BC_PERIODIC) NBC=FISHPAK_BC_PERIODIC

! Wannier flow boundaries

IF (PERIODIC_TEST==5) THEN
   LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
ENDIF

! Poisson solver with stretching in the 1st coordinate

SELECT_POISSON_SOLVER: SELECT CASE(M%IPS)

   CASE (0:1) SELECT_POISSON_SOLVER
      IF (.NOT.TWO_D) CALL H3CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,M%HX,XLM,ITRN,JTRN,IERR,M%SAVE1)
      IF (TWO_D .AND. .NOT.CYLINDRICAL) CALL H2CZIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HX,XLM,ITRN,IERR,M%SAVE1)
      IF (TWO_D .AND. CYLINDRICAL) THEN
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_DIRICHLET_DIRICHLET) LBC = 5
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_DIRICHLET_NEUMANN)   LBC = 6
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_NEUMANN_NEUMANN)     LBC = 6
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_NEUMANN_DIRICHLET)   LBC = 5
         CALL H2CYIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,XLM,XMU,ITRN,IERR,M%SAVE1)
      ENDIF
      IF (EXTERNAL_BOUNDARY_CORRECTION) THEN
         IF (.NOT.TWO_D) CALL H3CZIS(XS,XF,IBAR,LBC2,YS,YF,JBAR,MBC2,ZS,ZF,KBAR,NBC2,M%HX,XLM,ITRN,JTRN,IERR,M%SAVE2)
      ENDIF
   CASE (2) SELECT_POISSON_SOLVER
      CALL H3CZIS(YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HY,XLM,ITRN,JTRN,IERR,M%SAVE1)
   CASE (3) SELECT_POISSON_SOLVER
      IF (TWO_D) THEN
         CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE1)
      ELSE
         CALL H3CZIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,JTRN,IERR,M%SAVE1)
      ENDIF
   CASE (4) SELECT_POISSON_SOLVER
      CALL H3CSIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HX,M%HY)
   CASE (5) SELECT_POISSON_SOLVER
      IF (TWO_D) THEN
         CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE1)
      ELSE
         CALL H3CSIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HX,M%HZ)
      ENDIF
   CASE (6) SELECT_POISSON_SOLVER
      CALL H3CSIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HZ,M%HY)
   CASE (7) SELECT_POISSON_SOLVER
      CALL H2CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,M%HX,XLM,ITRN,IERR,M%SAVE1)

END SELECT SELECT_POISSON_SOLVER

! Specify the pressure boundary condition for each wall cell

WALL_CELL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS
   WC => M%WALL(IW)
   IOR = WC%ONE_D%IOR
   SELECT CASE(IOR)
      CASE( 1)
         IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. &
             LBC==FISHPAK_BC_NEUMANN_DIRICHLET   .OR. LBC==6) WC%PRESSURE_BC_INDEX = NEUMANN
         IF (LBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. &
             LBC==FISHPAK_BC_DIRICHLET_NEUMANN   .OR. LBC==5) WC%PRESSURE_BC_INDEX = DIRICHLET
      CASE(-1)
         IF (LBC==FISHPAK_BC_DIRICHLET_NEUMANN   .OR. &
             LBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. LBC==6) WC%PRESSURE_BC_INDEX = NEUMANN
         IF (LBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. &
             LBC==FISHPAK_BC_NEUMANN_DIRICHLET   .OR. LBC==5) WC%PRESSURE_BC_INDEX = DIRICHLET
      CASE( 2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. MBC==FISHPAK_BC_NEUMANN_DIRICHLET) WC%PRESSURE_BC_INDEX = NEUMANN
         IF (MBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. MBC==FISHPAK_BC_DIRICHLET_NEUMANN) WC%PRESSURE_BC_INDEX = DIRICHLET
      CASE(-2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. MBC==FISHPAK_BC_DIRICHLET_NEUMANN) WC%PRESSURE_BC_INDEX = NEUMANN
         IF (MBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. MBC==FISHPAK_BC_NEUMANN_DIRICHLET) WC%PRESSURE_BC_INDEX = DIRICHLET
      CASE( 3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. NBC==FISHPAK_BC_NEUMANN_DIRICHLET) WC%PRESSURE_BC_INDEX = NEUMANN
         IF (NBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. NBC==FISHPAK_BC_DIRICHLET_NEUMANN) WC%PRESSURE_BC_INDEX = DIRICHLET
      CASE(-3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. NBC==FISHPAK_BC_DIRICHLET_NEUMANN) WC%PRESSURE_BC_INDEX = NEUMANN
         IF (NBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. NBC==FISHPAK_BC_NEUMANN_DIRICHLET) WC%PRESSURE_BC_INDEX = DIRICHLET
   END SELECT
ENDDO WALL_CELL_LOOP

IF (EXTERNAL_BOUNDARY_CORRECTION) THEN
   WALL_CELL_LOOP_2: DO IW=1,M%N_EXTERNAL_WALL_CELLS
      WC => M%WALL(IW)
      IOR = WC%ONE_D%IOR
      WC%LAPLACE_BC_INDEX = NEUMANN
      SELECT CASE(IOR)
         CASE( 1)
            IF (LBC2==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. &
                LBC2==FISHPAK_BC_DIRICHLET_NEUMANN   .OR. LBC2==5) WC%LAPLACE_BC_INDEX = DIRICHLET
         CASE(-1)
            IF (LBC2==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. &
                LBC2==FISHPAK_BC_NEUMANN_DIRICHLET   .OR. LBC2==5) WC%LAPLACE_BC_INDEX = DIRICHLET
         CASE( 2)
            IF (MBC2==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. MBC2==FISHPAK_BC_DIRICHLET_NEUMANN) WC%LAPLACE_BC_INDEX = DIRICHLET
         CASE(-2)
            IF (MBC2==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. MBC2==FISHPAK_BC_NEUMANN_DIRICHLET) WC%LAPLACE_BC_INDEX = DIRICHLET
         CASE( 3)
            IF (NBC2==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. NBC2==FISHPAK_BC_DIRICHLET_NEUMANN) WC%LAPLACE_BC_INDEX = DIRICHLET
         CASE(-3)
            IF (NBC2==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. NBC2==FISHPAK_BC_NEUMANN_DIRICHLET) WC%LAPLACE_BC_INDEX = DIRICHLET
      END SELECT
   ENDDO WALL_CELL_LOOP_2
ENDIF

! Check for errors with Poisson solver initialization

IF (IERR/=0) THEN
   WRITE(LU_ERR,'(A,I2,A,I3)') 'ERROR: Poisson initialization error, Number=',IERR, ', Mesh=',NM
   STOP_STATUS = SETUP_STOP
   RETURN
ENDIF

END SUBROUTINE INITIALIZE_POISSON_SOLVER


!SUBROUTINE INITIALIZE_LAPLACE_SOLVER
!USE PRES, ONLY: BUILD_SPARSE_MATRIX_LAPLACE
!USE MKL_DSS
!TYPE(MKL_DSS_HANDLE) :: handle
!INTEGER :: IERR,perm(1)
!
!IF (.NOT.LAPLACE_PRESSURE_CORRECTION) RETURN
!
!CALL BUILD_SPARSE_MATRIX_LAPLACE(NM)
!
!! Initialize the solver
!IERR = dss_create( handle, MKL_DSS_DEFAULTS )
!
!! Define the non-zero structure of the matrix
!IERR = dss_define_structure( handle, MKL_DSS_SYMMETRIC, M%A_ROW_INDEX, M%A_N_ROWS, M%A_N_COLS, M%A_COLUMNS, M%A_N_ELEMENTS )
!print *,'dss_define_structure complete...'
!
!! Reorder the matrix
!perm(1)=0
!IERR = dss_reorder( handle, MKL_DSS_DEFAULTS, perm )
!print *,'dss_reorder complete...'
!
!! Factor the matrix
!IERR = dss_factor_real( handle, MKL_DSS_DEFAULTS, M%A_VALUES )
!print *,'dss_factor_real complete...'
!
!stop
!
!END SUBROUTINE INITIALIZE_LAPLACE_SOLVER


SUBROUTINE INITIALIZE_INTERPOLATION

! Create arrays by which info is to exchanged across meshes

INTEGER :: NOM,I,J,K
TYPE (MESH_TYPE), POINTER :: M2

ALLOCATE(M%INTERPOLATED_MESH(1:M%IBAR,1:M%JBAR,1:M%KBAR),  STAT=IZERO)
CALL ChkMemErr('INIT','INTERPOLATED_MESH',IZERO)
M%INTERPOLATED_MESH = 0

IF (EVACUATION_ONLY(NM)) RETURN

DO K=1,M%KBAR
   DO J=1,M%JBAR
      DO I=1,M%IBAR
         OTHER_MESH_LOOP: DO NOM=1,NM-1
            M2=>MESHES(NOM)
            IF (M%X(I-1)>=M2%XS .AND. M%X(I)<=M2%XF .AND. &
                M%Y(J-1)>=M2%YS .AND. M%Y(J)<=M2%YF .AND. &
                M%Z(K-1)>=M2%ZS .AND. M%Z(K)<=M2%ZF) THEN
               M%INTERPOLATED_MESH(I,J,K) = NOM
               EXIT OTHER_MESH_LOOP
            ENDIF
         ENDDO OTHER_MESH_LOOP
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE INITIALIZE_INTERPOLATION

END SUBROUTINE INITIALIZE_MESH_VARIABLES_2



SUBROUTINE INITIALIZE_DEVICES(NM)

USE COMPLEX_GEOMETRY, ONLY : GET_CFACE_INDEX

! Find the WALL_INDEX for a device that is near a solid wall

INTEGER, INTENT(IN) :: NM
INTEGER :: III,N,II,JJ,KK,IOR,IW,SURF_INDEX,IIG,JJG,KKG,ICF
REAL(EB) :: DEPTH
TYPE (DEVICE_TYPE), POINTER :: DV
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)

DEVICE_LOOP: DO N=1,N_DEVC

   DV => DEVICE(N)

   IF (DV%QUANTITY_INDEX(1)>=0) CYCLE DEVICE_LOOP  ! Do not process gas phsae devices

   IF (DV%INIT_ID=='null') THEN ! Assume the device is tied to a WALL cell or CFACE

      IF (NM/=DV%MESH) CYCLE DEVICE_LOOP
      II  = INT(GINV(DV%X-M%XS,1,NM)*M%RDXI   + 1._EB)
      JJ  = INT(GINV(DV%Y-M%YS,2,NM)*M%RDETA  + 1._EB)
      KK  = INT(GINV(DV%Z-M%ZS,3,NM)*M%RDZETA + 1._EB)
      IIG = II
      JJG = JJ
      KKG = KK
      IOR = DV%IOR
      IW  = 0
      ICF = 0

      IF (TRIM(DV%QUANTITY(1))=='SOLID CELL TEMPERATURE') THEN
         ! For SOLID CELL TEMPERATURE (II,JJ,KK) should be inside SOLID,
         ! our task is to find the first gas phase cell (IIG,JJG,KKG) in direction IOR,
         ! currently assumes II and IIG, etc., are on the same mesh
         SELECT CASE (IOR)
            CASE ( 1); DO IIG=II,IBP1; IF (.NOT.M%SOLID(M%CELL_INDEX(IIG,JJG,KKG))) EXIT; ENDDO
            CASE (-1); DO IIG=II,0,-1; IF (.NOT.M%SOLID(M%CELL_INDEX(IIG,JJG,KKG))) EXIT; ENDDO
            CASE ( 2); DO JJG=JJ,JBP1; IF (.NOT.M%SOLID(M%CELL_INDEX(IIG,JJG,KKG))) EXIT; ENDDO
            CASE (-2); DO JJG=JJ,0,-1; IF (.NOT.M%SOLID(M%CELL_INDEX(IIG,JJG,KKG))) EXIT; ENDDO
            CASE ( 3); DO KKG=KK,KBP1; IF (.NOT.M%SOLID(M%CELL_INDEX(IIG,JJG,KKG))) EXIT; ENDDO
            CASE (-3); DO KKG=KK,0,-1; IF (.NOT.M%SOLID(M%CELL_INDEX(IIG,JJG,KKG))) EXIT; ENDDO
         END SELECT
      ENDIF

      IF (IOR/=0) CALL GET_WALL_INDEX(NM,IIG,JJG,KKG,IOR,IW)

      IF (IW==0)  CALL GET_CFACE_INDEX(NM,IIG,JJG,KKG,DV%X,DV%Y,DV%Z,ICF)

      IF (IW==0 .AND. ICF==0 .AND. DV%SPATIAL_STATISTIC=='null') THEN
         WRITE(LU_ERR,'(A,I0,A,A)') 'ERROR: Reposition DEVC No.',DV%ORDINAL,', ID = ',TRIM(DV%ID)
         WRITE(LU_ERR,'(A)') 'FDS cannot determine which boundary cell to assign.'
         STOP_STATUS = SETUP_STOP
         RETURN
      ELSEIF (IW>0) THEN
         DV%WALL_INDEX = IW
         SURF_INDEX = M%WALL(IW)%SURF_INDEX
      ELSEIF (ICF>0) THEN
         DV%CFACE_INDEX = ICF
         SURF_INDEX = M%CFACE(ICF)%SURF_INDEX
      ELSE
         SURF_INDEX = DV%SURF_INDEX
      ENDIF

   ELSE ! Assume the device is tied to a particle

      IF (DV%PART_CLASS_INDEX<1) CYCLE DEVICE_LOOP
      SURF_INDEX = LAGRANGIAN_PARTICLE_CLASS(DV%PART_CLASS_INDEX)%SURF_INDEX

   ENDIF

   ! Make sure that thermally-thick output is appropriate

   IF (OUTPUT_QUANTITY(DV%QUANTITY_INDEX(1))%INSIDE_SOLID) THEN
      IF (SURFACE(SURF_INDEX)%THERMAL_BC_INDEX /= THERMALLY_THICK) THEN
         WRITE(LU_ERR,'(A,I3,A)') 'ERROR: DEVC ',N,' must be associated with a heat-conducting surface'
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF
      IF (DV%DEPTH>TWO_EPSILON_EB) THEN
         DEPTH = DV%DEPTH
      ELSE
         DEPTH = MAX(0._EB,SUM(SURFACE(SURF_INDEX)%LAYER_THICKNESS)+DV%DEPTH)
      ENDIF
      DV%I_DEPTH = SURFACE(SURF_INDEX)%N_CELLS_INI
      DO III=SURFACE(SURF_INDEX)%N_CELLS_INI,1,-1
         IF (DEPTH<=SURFACE(SURF_INDEX)%X_S(III)) DV%I_DEPTH = III
      ENDDO
   ENDIF

ENDDO DEVICE_LOOP

END SUBROUTINE INITIALIZE_DEVICES


SUBROUTINE INITIALIZE_PROFILES(NM)

INTEGER, INTENT(IN) :: NM
INTEGER :: NN,N,II,JJ,KK,IW,IOR
LOGICAL :: SUCCESS
TYPE (PROFILE_TYPE), POINTER :: PF
TYPE (MESH_TYPE), POINTER :: M
TYPE (SURFACE_TYPE), POINTER :: SF

IF (EVACUATION_ONLY(NM)) RETURN

M => MESHES(NM)

PROF_LOOP: DO N=1,N_PROF

   PF => PROFILE(N)
   IF (NM/=PF%MESH) CYCLE PROF_LOOP
   II  = INT(GINV(PF%X-M%XS,1,NM)*M%RDXI   + 1._EB)
   JJ  = INT(GINV(PF%Y-M%YS,2,NM)*M%RDETA  + 1._EB)
   KK  = INT(GINV(PF%Z-M%ZS,3,NM)*M%RDZETA + 1._EB)

   IF (PF%IOR/=0) THEN

      ! The PROFile is for a WALL cell

      IOR = PF%IOR
      CALL GET_WALL_INDEX(NM,II,JJ,KK,IOR,IW)
      IF (IW>0) THEN
         PF%WALL_INDEX = IW
         SF => SURFACE(M%WALL(IW)%SURF_INDEX)
      ELSE
         WRITE(LU_ERR,'(A,I0,A)') 'ERROR: Reposition PROF No. ',PF%ORDINAL, '. FDS cannot determine which boundary cell to assign'
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF

   ELSE  ! The PROFile is for a Lagrangian PARTicle

      SF => SURFACE(LAGRANGIAN_PARTICLE_CLASS(PF%PART_CLASS_INDEX)%SURF_INDEX)

   ENDIF

   IF (SF%THERMAL_BC_INDEX/=THERMALLY_THICK) THEN
      WRITE(LU_ERR,'(A,I3,A)') 'ERROR: PROFile ',N, ' must be associated with a heat-conducting surface'
      STOP_STATUS = SETUP_STOP
      RETURN
   ENDIF

   IF (PF%QUANTITY /= 'TEMPERATURE' .AND. PF%QUANTITY /= 'DENSITY') THEN
      SUCCESS = .FALSE.
      DO NN=1,SF%N_MATL
         IF (PF%QUANTITY==SF%MATL_NAME(NN)) SUCCESS = .TRUE.
      ENDDO
      IF (.NOT.SUCCESS) THEN
         WRITE(LU_ERR,'(A,A,A)') 'ERROR: QUANTITY ',TRIM(PF%QUANTITY), ' is not appropriate for the designated location'
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF
   ENDIF

ENDDO PROF_LOOP

END SUBROUTINE INITIALIZE_PROFILES


SUBROUTINE GET_WALL_INDEX(NM,II,JJ,KK,IOR,IW)

! Find the wall index corresponding to the -IOR face of cell (II,JJ,KK)

INTEGER, INTENT(IN) :: NM,IOR
INTEGER, INTENT(OUT) :: IW
INTEGER :: IC,II,JJ,KK
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
IC  = M%CELL_INDEX(II,JJ,KK)

IF (M%SOLID(IC)) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      II = II-1
      CASE( 1)
         IF (II<M%IBP1) II = II+1
      CASE(-2)
         IF (JJ>0)      JJ = JJ-1
      CASE( 2)
         IF (JJ<M%JBP1) JJ = JJ+1
      CASE(-3)
         IF (KK>0)      KK = KK-1
      CASE( 3)
         IF (KK<M%KBP1) KK = KK+1
   END SELECT
ENDIF

IC  = M%CELL_INDEX(II,JJ,KK)
IW  = M%WALL_INDEX(IC,-IOR)

IF (IW<=0) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      IC = M%CELL_INDEX(II-1,JJ,KK)
      CASE( 1)
         IF (II<M%IBP1) IC = M%CELL_INDEX(II+1,JJ,KK)
      CASE(-2)
         IF (JJ>0)      IC = M%CELL_INDEX(II,JJ-1,KK)
      CASE( 2)
         IF (JJ<M%JBP1) IC = M%CELL_INDEX(II,JJ+1,KK)
      CASE(-3)
         IF (KK>0)      IC = M%CELL_INDEX(II,JJ,KK-1)
      CASE( 3)
         IF (KK<M%KBP1) IC = M%CELL_INDEX(II,JJ,KK+1)
   END SELECT
   IW = M%WALL_INDEX(IC,-IOR)
ENDIF

END SUBROUTINE GET_WALL_INDEX



SUBROUTINE INITIALIZE_GLOBAL_VARIABLES

USE CONTROL_VARIABLES, ONLY: N_CTRL
INTEGER :: IZERO, IG

! Initialize time, printout and plot clocks

ALLOCATE(PART_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PART_CLOCK',IZERO)
ALLOCATE(ISOF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','ISOF_CLOCK',IZERO)
ALLOCATE(BNDF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','BNDF_CLOCK',IZERO)
ALLOCATE(SLCF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','SLCF_CLOCK',IZERO)
ALLOCATE(RADF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','RADF_CLOCK',IZERO)
ALLOCATE(SL3D_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','SLCF_CLOCK',IZERO)
ALLOCATE(PL3D_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PL3D_CLOCK',IZERO)
ALLOCATE(PROF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PROF_CLOCK',IZERO)

ICYC              = 0
PART_CLOCK        = T_BEGIN
DEVC_CLOCK        = T_BEGIN
CTRL_CLOCK        = T_BEGIN
PROF_CLOCK        = T_BEGIN
PL3D_CLOCK        = T_BEGIN + DT_PL3D
ISOF_CLOCK        = T_BEGIN
SLCF_CLOCK        = T_BEGIN
RADF_CLOCK        = T_RADF_BEGIN
SL3D_CLOCK        = T_BEGIN + DT_SL3D
GEOC_CLOCK        = T_BEGIN
GEOM_CLOCK        = T_BEGIN
BNDC_CLOCK        = T_BEGIN
BNDF_CLOCK        = T_BEGIN
RESTART_CLOCK     = T_BEGIN + DT_RESTART
FLUSH_CLOCK       = T_BEGIN + DT_FLUSH
HRR_CLOCK         = T_BEGIN
CPU_CLOCK         = T_BEGIN + DT_CPU
MINT_CLOCK        = T_BEGIN
T_LAST_DUMP_HRR   = T_BEGIN
T_LAST_DUMP_MASS  = T_BEGIN
T_LAST_DUMP_MOM   = T_BEGIN

IF (N_DEVC==0)      DEVC_CLOCK = 1.E10_EB
IF (N_CTRL==0)      CTRL_CLOCK = 1.E10_EB
IF (N_PROF==0)      PROF_CLOCK = 1.E10_EB
IF (N_ISOF==0)      ISOF_CLOCK = 1.E10_EB
! N_FACE manages the geometry output time GEOM_CLOCK:
DO IG=1,N_GEOMETRY; N_FACE = N_FACE + GEOMETRY(IG)%N_FACES; ENDDO
IF (N_FACE==0)      GEOC_CLOCK = 1.E10_EB
IF (N_FACE==0)      GEOM_CLOCK = 1.E10_EB
IF (N_FACE==0)      BNDC_CLOCK = 1.E10_EB
IF (N_BNDF==0)      BNDF_CLOCK = 1.E10_EB
IF (.NOT.MASS_FILE) MINT_CLOCK = 1.E10_EB

ALLOCATE(ENTHALPY_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','ENTHALPY_SUM',IZERO)
ENTHALPY_SUM = 0._EB
ALLOCATE(Q_DOT(N_Q_DOT,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','Q_DOT',IZERO)
Q_DOT = 0._EB
ALLOCATE(Q_DOT_SUM(N_Q_DOT,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','Q_DOT_SUM',IZERO)
Q_DOT_SUM = 0._EB
ALLOCATE(M_DOT(N_M_DOT,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','M_DOT',IZERO)
M_DOT = 0._EB
ALLOCATE(M_DOT_SUM(N_M_DOT,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','M_DOT_SUM',IZERO)
M_DOT_SUM=0._EB

ALLOCATE(MINT(0:N_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT',IZERO)
ALLOCATE(MINT_SUM(0:N_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT_SUM',IZERO)
MINT_SUM=0._EB

ALLOCATE(PRESSURE_ERROR_MAX(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_ERROR_MAX',IZERO)
ALLOCATE(PRESSURE_ERROR_MAX_LOC(3,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_ERROR_MAX_LOC',IZERO)
PRESSURE_ERROR_MAX     = 0._EB
PRESSURE_ERROR_MAX_LOC = 0

ALLOCATE(VELOCITY_ERROR_MAX(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX',IZERO)
ALLOCATE(VELOCITY_ERROR_MAX_LOC(3,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX_LOC',IZERO)
VELOCITY_ERROR_MAX     = 0._EB
VELOCITY_ERROR_MAX_LOC = 0

END SUBROUTINE INITIALIZE_GLOBAL_VARIABLES


SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)

! Initialize wall cell variables at external and obstruction boundaries

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE
USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT
INTEGER, INTENT(IN) :: I,J,K,NM,OBST_INDEX,IW,IOR,SURF_INDEX
INTEGER  :: NOM_FOUND,NOM=0,ITER,IIO_MIN,IIO_MAX,JJO_MIN,JJO_MAX,KKO_MIN,KKO_MAX,VENT_INDEX
INTEGER, INTENT(OUT) :: IERR
REAL(EB), INTENT(IN) :: TT
REAL(EB) :: PX,PY,PZ,T_ACTIVATE,XIN,YIN,ZIN,DIST,XW,YW,ZW,UW,RDN,AW,TSI,&
            ZZ_GET(1:N_TRACKED_SPECIES),RSUM_F,R1,RR,DELTA
INTEGER  :: N,SURF_INDEX_NEW,IIG,JJG,KKG,IIO,JJO,KKO,IC,ICG,ICO,NOM_CHECK(0:1),BOUNDARY_TYPE
INTEGER :: NSLICE
LOGICAL :: VENT_FOUND,ALIGNED
TYPE (MESH_TYPE), POINTER :: M,MM=>NULL()
TYPE (OBSTRUCTION_TYPE), POINTER :: OBX=>NULL()
TYPE (VENTS_TYPE), POINTER :: VT=>NULL()
TYPE (WALL_TYPE), POINTER :: WC=>NULL()
TYPE (SURFACE_TYPE), POINTER :: SF

IERR = 0
M=>MESHES(NM)

! Determine if a VENT covers the surface

VENT_INDEX = 0
SURF_INDEX_NEW = SURF_INDEX
VENT_FOUND = .FALSE.

VENT_SEARCH_LOOP: DO N=1,M%N_VENT

   VT => M%VENTS(N)
   IF (OBST_INDEX>0) THEN
      IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY)           CYCLE VENT_SEARCH_LOOP
      IF (.NOT.M%OBSTRUCTION(OBST_INDEX)%ALLOW_VENT) CYCLE VENT_SEARCH_LOOP
      IF (VT%OBST_INDEX>0 .AND. VT%OBST_INDEX/=OBST_INDEX) CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (VT%IOR/=IOR) CYCLE VENT_SEARCH_LOOP

   IF (ABS(IOR)==1) THEN
      IF (IOR== 1 .AND. I/=VT%I1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-1 .AND. I/=VT%I1+1) CYCLE VENT_SEARCH_LOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VENT_SEARCH_LOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VENT_SEARCH_LOOP
      IF ( VT%RADIUS>0._EB .AND. ((M%YC(J)-VT%Y0)**2 + (M%ZC(K)-VT%Z0)**2)>(VT%RADIUS**2) ) CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (ABS(IOR)==2) THEN
      IF (IOR== 2 .AND. J/=VT%J1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-2 .AND. J/=VT%J1+1) CYCLE VENT_SEARCH_LOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VENT_SEARCH_LOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VENT_SEARCH_LOOP
      IF ( VT%RADIUS>0._EB .AND. ((M%XC(I)-VT%X0)**2 + (M%ZC(K)-VT%Z0)**2)>(VT%RADIUS**2) ) CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (ABS(IOR)==3) THEN
      IF (IOR== 3 .AND. K/=VT%K1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-3 .AND. K/=VT%K1+1) CYCLE VENT_SEARCH_LOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VENT_SEARCH_LOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VENT_SEARCH_LOOP
      IF ( VT%RADIUS>0._EB .AND. ((M%XC(I)-VT%X0)**2 + (M%YC(J)-VT%Y0)**2)>(VT%RADIUS**2) ) CYCLE VENT_SEARCH_LOOP
   ENDIF

   ! Check if there are over-lapping VENTs

   IF (VENT_FOUND) THEN
      WRITE(LU_ERR,'(A,I0,A,3(I0,1X),A,I0,A)') 'WARNING: Two VENTs overlap in MESH ',NM,', Cell ',I,J,K,'. VENT ', &
                                           VT%ORDINAL,' rejected for that cell'
      EXIT VENT_SEARCH_LOOP
   ENDIF

   VENT_FOUND = .TRUE.

   ! Reassign the SURF index to be that of the VENT

   VENT_INDEX = N
   SURF_INDEX_NEW = VT%SURF_INDEX

ENDDO VENT_SEARCH_LOOP

! Compute boundary cell physical coords (XW,YW,ZW) and area (AW)

IF (ABS(IOR)==1) THEN
   IF (IOR== 1) THEN
      XW     = M%X(I)
      IIG    = I+1
      RDN    = M%RDXN(I)
      AW     = M%R(I)*M%DY(J)*M%DZ(K)
      UW     = -U0
   ENDIF
   IF (IOR==-1) THEN
      XW     = M%X(I-1)
      IIG    = I-1
      RDN    = M%RDXN(I-1)
      AW     = M%R(I-1)*M%DY(J)*M%DZ(K)
      UW     = U0
   ENDIF
   JJG    = J
   KKG    = K
   YW = M%YC(J)
   ZW = M%ZC(K)
ENDIF
IF (ABS(IOR)==2) THEN
   IF (IOR== 2) THEN
      YW     = M%Y(J)
      JJG    = J+1
      RDN    = M%RDYN(J)
      UW     = -V0
   ENDIF
   IF (IOR==-2) THEN
      YW     = M%Y(J-1)
      JJG    = J-1
      RDN    = M%RDYN(J-1)
      UW     = V0
   ENDIF
   IIG    = I
   KKG    = K
   XW = M%XC(I)
   ZW = M%ZC(K)
   AW = M%DX(I)*M%DZ(K)
ENDIF
IF (ABS(IOR)==3) THEN
   IF (IOR== 3) THEN
      ZW     = M%Z(K)
      KKG    = K+1
      RDN    = M%RDZN(K)
      UW     = -W0
   ENDIF
   IF (IOR==-3) THEN
      ZW     = M%Z(K-1)
      KKG    = K-1
      RDN    = M%RDZN(K-1)
      UW     = W0
   ENDIF
   IIG    = I
   JJG    = J
   XW = M%XC(I)
   YW = M%YC(J)
   AW = M%DX(I)*M%RC(I)*M%DY(J)
ENDIF

IF (IOR==0) THEN
   IIG = I
   JJG = J
   KKG = K
ENDIF

! Save the wall index

IC  = M%CELL_INDEX(I  ,J  ,K  )
ICG = M%CELL_INDEX(IIG,JJG,KKG)
M%WALL_INDEX(ICG,-IOR) = IW

! Use BOUNDARY_TYPE to indicate whether the boundary cell is blocked or on an obstruction that is HIDDEN

BOUNDARY_TYPE = NULL_BOUNDARY

IF (IW<=M%N_EXTERNAL_WALL_CELLS .AND. OBST_INDEX==0) BOUNDARY_TYPE = SOLID_BOUNDARY

IF (OBST_INDEX>0) THEN
   IF (.NOT.M%OBSTRUCTION(OBST_INDEX)%HIDDEN) THEN
      BOUNDARY_TYPE = SOLID_BOUNDARY
      IF (IW<=M%N_EXTERNAL_WALL_CELLS) M%SOLID(IC) = .TRUE.
   ENDIF
ENDIF

IF (M%SOLID(ICG)) BOUNDARY_TYPE = NULL_BOUNDARY

! Check for neighboring meshes in a multiple mesh calculation

NOM_FOUND = 0
IIO_MIN =  1000000
IIO_MAX = -1000000
JJO_MIN =  1000000
JJO_MAX = -1000000
KKO_MIN =  1000000
KKO_MAX = -1000000
NOM_CHECK = 0

CHECK_MESHES: IF (IW<=M%N_EXTERNAL_WALL_CELLS .AND. .NOT.EVACUATION_ONLY(NM)) THEN

   DO ITER=0,1
      XIN = XW
      YIN = YW
      ZIN = ZW
      IF (SURF_INDEX_NEW==PERIODIC_SURF_INDEX .OR. SURF_INDEX_NEW==PERIODIC_FLOW_ONLY_SURF_INDEX) THEN
         SELECT CASE(IOR)
            CASE( 1) ; XIN = XF_MAX
            CASE(-1) ; XIN = XS_MIN
            CASE( 2) ; YIN = YF_MAX
            CASE(-2) ; YIN = YS_MIN
            CASE( 3) ; ZIN = ZF_MAX
            CASE(-3) ; ZIN = ZS_MIN
         END SELECT
      ENDIF
      IF (ABS(IOR)/=1) XIN = XW + (ITER*0.95_EB-0.475_EB)*(M%X(I)-M%X(I-1))
      IF (ABS(IOR)/=2) YIN = YW + (ITER*0.95_EB-0.475_EB)*(M%Y(J)-M%Y(J-1))
      IF (ABS(IOR)/=3) ZIN = ZW + (ITER*0.95_EB-0.475_EB)*(M%Z(K)-M%Z(K-1))
      IF (IOR== 1) XIN = XIN - MESH_SEPARATION_DISTANCE
      IF (IOR==-1) XIN = XIN + MESH_SEPARATION_DISTANCE
      IF (IOR== 2) YIN = YIN - MESH_SEPARATION_DISTANCE
      IF (IOR==-2) YIN = YIN + MESH_SEPARATION_DISTANCE
      IF (IOR== 3) ZIN = ZIN - MESH_SEPARATION_DISTANCE
      IF (IOR==-3) ZIN = ZIN + MESH_SEPARATION_DISTANCE
      CALL SEARCH_OTHER_MESHES(XIN,YIN,ZIN,NOM,IIO,JJO,KKO)
      NOM_CHECK(ITER) = NOM
      IF (NOM/=0) THEN
         IIO_MIN = MIN(IIO_MIN,IIO)
         IIO_MAX = MAX(IIO_MAX,IIO)
         JJO_MIN = MIN(JJO_MIN,JJO)
         JJO_MAX = MAX(JJO_MAX,JJO)
         KKO_MIN = MIN(KKO_MIN,KKO)
         KKO_MAX = MAX(KKO_MAX,KKO)
      ENDIF
   ENDDO

   ! Check to see if the current interpolated cell face spans more than one other mesh

   IF (NOM_CHECK(0)/=NOM_CHECK(1)) THEN
      WRITE(LU_ERR,'(A,I3,A,I3)') 'ERROR: MESH ',NM,' is not in alignment with MESH ',MAXVAL(NOM_CHECK)
      STOP_STATUS = SETUP_STOP
      IERR = 1
      RETURN
   ENDIF

   FOUND_OTHER_MESH: IF (NOM>0) THEN
      MM=>MESHES(NOM)
      ALIGNED = .TRUE.
      IF ( (ABS(IOR)==2 .OR. ABS(IOR)==3) .AND. MM%DX(IIO_MIN)<=M%DX(I) ) THEN
         IF (ABS( ((MM%X(IIO_MAX)-MM%X(IIO_MIN-1))-(M%X(I)-M%X(I-1))) / MM%DX(IIO_MIN))>ALIGNMENT_TOLERANCE ) ALIGNED = .FALSE.
      ENDIF
      IF ( (ABS(IOR)==1 .OR. ABS(IOR)==3) .AND. MM%DY(JJO_MIN)<=M%DY(J) ) THEN
         IF (ABS( ((MM%Y(JJO_MAX)-MM%Y(JJO_MIN-1))-(M%Y(J)-M%Y(J-1))) / MM%DY(JJO_MIN))>ALIGNMENT_TOLERANCE ) ALIGNED = .FALSE.
      ENDIF
      IF ( (ABS(IOR)==1 .OR. ABS(IOR)==2) .AND. MM%DZ(KKO_MIN)<=M%DZ(K) ) THEN
         IF (ABS( ((MM%Z(KKO_MAX)-MM%Z(KKO_MIN-1))-(M%Z(K)-M%Z(K-1))) / MM%DZ(KKO_MIN))>ALIGNMENT_TOLERANCE ) ALIGNED = .FALSE.
      ENDIF
      IF (.NOT.ALIGNED) THEN
         WRITE(LU_ERR,'(A,I3,A,I3)') 'ERROR: MESH ',NM,' is out of alignment with MESH ',NOM
         STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF

      SELECT CASE(ABS(IOR))
         CASE(1)
            IF ( (M%DY(J)<0.99_EB*MM%DY(JJO_MIN)) .AND. (M%DZ(K)>1.01_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DY(J)>1.01_EB*MM%DY(JJO_MIN)) .AND. (M%DZ(K)<0.99_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
         CASE(2)
            IF ( (M%DX(I)<0.99_EB*MM%DX(IIO_MIN)) .AND. (M%DZ(K)>1.01_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DX(I)>1.01_EB*MM%DX(IIO_MIN)) .AND. (M%DZ(K)<0.99_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
         CASE(3)
            IF ( (M%DY(J)<0.99_EB*MM%DY(JJO_MIN)) .AND. (M%DX(I)>1.01_EB*MM%DX(IIO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DY(J)>1.01_EB*MM%DY(JJO_MIN)) .AND. (M%DX(I)<0.99_EB*MM%DX(IIO_MIN)) ) ALIGNED = .FALSE.
      END SELECT
      IF (.NOT.ALIGNED) THEN
         WRITE(LU_ERR,'(A,I3,A,I3)') 'ERROR: MESH ',NM,' is finer in one direction and coarser in the other than MESH ',NOM
         STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF

      ! NOM = "Number of the Other Mesh"
      !
      !     Mesh 1       Mesh 2
      !      3x6          1x2
      !  -------------------------
      !  |   |   |   |           |
      !  |---|---|---|           |
      !  |   |   |   #->         |
      !  |---|---|---|           |
      !  |   |   | <-*           |
      !  |---|---|---|-----------|
      !  |   |   |   |           |
      !  |---|---|---|           |
      !  |   |   |   |           |
      !  |---|---|---|           |
      !  |   |   |   |           |
      !  -------------------------
      !
      !                         NM=1,NOM=2,IW=*        NM=2,NOM=1,IW=#

      NOM_FOUND = NOM
      ICO = MM%CELL_INDEX(IIO_MIN,JJO_MIN,KKO_MIN)

      IF (OBST_INDEX==0) THEN
         IF (.NOT.M%SOLID(ICG) .AND. .NOT.MM%SOLID(ICO)) THEN
            BOUNDARY_TYPE = INTERPOLATED_BOUNDARY
            IF (SURF_INDEX_NEW/=PERIODIC_FLOW_ONLY_SURF_INDEX) SURF_INDEX_NEW = INTERPOLATED_SURF_INDEX
         ENDIF
         IF (M%SOLID(ICG) .OR. MM%SOLID(ICO)) THEN
            IF (MM%SOLID(ICO)) THEN
               SURF_INDEX_NEW = MM%OBSTRUCTION(MM%OBST_INDEX_C(ICO))%SURF_INDEX(IOR)
            ELSE
               SURF_INDEX_NEW = M%OBSTRUCTION(M%OBST_INDEX_C(ICG))%SURF_INDEX(-IOR)
            ENDIF
         ENDIF
         IF (M%SOLID(ICG) .AND. MM%SOLID(ICO) .AND.  &
            M%OBSTRUCTION(M%OBST_INDEX_C(ICG))%HT3D .AND.  MM%OBSTRUCTION(MM%OBST_INDEX_C(ICO))%HT3D) THEN
            SURF_INDEX_NEW = INTERPOLATED_SURF_INDEX
            BOUNDARY_TYPE = NULL_BOUNDARY
         ENDIF
         VENT_INDEX = 0
      ENDIF

      ! Do not allow a MIRROR boundary to sit on a mesh interface

      IF (VENT_INDEX>0) THEN
         IF (M%VENTS(VENT_INDEX)%BOUNDARY_TYPE==MIRROR_BOUNDARY) VENT_INDEX = 0
      ENDIF

      ! Open up the ghost cell at the interpolated boundary

      IF (BOUNDARY_TYPE == INTERPOLATED_BOUNDARY) M%SOLID(M%CELL_INDEX(I,J,K)) = .FALSE.

   ENDIF FOUND_OTHER_MESH

ENDIF CHECK_MESHES

! Initialize solid properties included in the WALL derived type variable

CALL ALLOCATE_STORAGE(NM,SURF_INDEX_NEW,WALL_INDEX=IW)

! Initialize wall cell (WC) variables

SF => SURFACE(SURF_INDEX_NEW)
WC => M%WALL(IW)

WC%ONE_D%II = I
WC%ONE_D%JJ = J
WC%ONE_D%KK = K
WC%ONE_D%IOR = IOR
WC%SURF_INDEX = SURF_INDEX_NEW
WC%OBST_INDEX = OBST_INDEX
WC%BOUNDARY_TYPE = BOUNDARY_TYPE

WC%X = XW
WC%Y = YW
WC%Z = ZW
WC%ONE_D%IIG    = IIG
WC%ONE_D%JJG    = JJG
WC%ONE_D%KKG    = KKG
WC%ONE_D%RDN    = RDN
WC%ONE_D%AREA   = AW
WC%ONE_D%U_NORMAL   = UW
WC%ONE_D%U_NORMAL_S = UW

IF (IW<=M%N_EXTERNAL_WALL_CELLS) THEN
   M%EXTERNAL_WALL(IW)%NOM     = NOM_FOUND
   M%EXTERNAL_WALL(IW)%IIO_MIN = IIO_MIN
   M%EXTERNAL_WALL(IW)%JJO_MIN = JJO_MIN
   M%EXTERNAL_WALL(IW)%KKO_MIN = KKO_MIN
   M%EXTERNAL_WALL(IW)%IIO_MAX = IIO_MAX
   M%EXTERNAL_WALL(IW)%JJO_MAX = JJO_MAX
   M%EXTERNAL_WALL(IW)%KKO_MAX = KKO_MAX
ENDIF

! If the simulation is only a TGA analysis, get the wall index

IF (WC%SURF_INDEX==TGA_SURF_INDEX) TGA_WALL_INDEX = IW

! Fill array containing K index of terrain following slice, used in dump slice

DO NSLICE = 1, M%N_TERRAIN_SLCF
   IF (IOR==3) M%K_AGL_SLICE(I,J,NSLICE) = M%K_AGL_SLICE(I,J,NSLICE) + WC%ONE_D%KKG
ENDDO

! Do not assign normal velocities at boundaries of evacuation meshes

IF (EVACUATION_ONLY(NM)) WC%ONE_D%U_NORMAL = 0._EB

! Assign internal values of temp, density, and mass fraction

WC%ONE_D%RHO_F = M%RHO(IIG,JJG,KKG)
WC%ONE_D%U_TAU = 0._EB
WC%ONE_D%Y_PLUS = 1._EB
WC%ONE_D%Z_STAR = 1._EB
WC%ONE_D%RHO_D_F = 0._EB
WC%ONE_D%RHO_D_DZDN_F = 0._EB

IF (.NOT.EVACUATION_ONLY(NM) .AND. (M%SOLID(M%CELL_INDEX(I,J,K)).OR.IW<=M%N_EXTERNAL_WALL_CELLS)) THEN
   M%RSUM(I,J,K) = M%RSUM(IIG,JJG,KKG)
   WC%ONE_D%ZZ_F(1:N_TRACKED_SPECIES)  = M%ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES)
   M%ZZ(I,J,K,1:N_TRACKED_SPECIES) = M%ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES)
ENDIF

! Compute the mass of the grid cell corresponding to the wall cell

IF (OBST_INDEX>0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   OBX=>M%OBSTRUCTION(OBST_INDEX)
   WC%ONE_D%AREA_ADJUST = SF%AREA_MULTIPLIER*OBX%INPUT_AREA(ABS(IOR))/OBX%FDS_AREA(ABS(IOR))
   IF (WC%ONE_D%AREA_ADJUST<=TWO_EPSILON_EB) WC%ONE_D%AREA_ADJUST = 1._EB
   IF (OBX%MASS>1.E5_EB) THEN
      OBX%MASS = SF%SURFACE_DENSITY*WC%ONE_D%AREA*WC%ONE_D%AREA_ADJUST
      IF (OBX%BULK_DENSITY>0._EB) OBX%MASS = OBX%BULK_DENSITY*(OBX%X2-OBX%X1)*(OBX%Y2-OBX%Y1)*(OBX%Z2-OBX%Z1)/OBX%VOLUME_ADJUST
   ENDIF
ENDIF

! Prescribe normal velocity for surface cell (U_NORMAL_0)

WC%ONE_D%U_NORMAL_0 = SF%VEL

IF (OBST_INDEX>0 .AND. ABS(SF%VOLUME_FLOW)>=TWO_EPSILON_EB .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   OBX=>M%OBSTRUCTION(OBST_INDEX)
   WC%ONE_D%U_NORMAL_0 = SF%VOLUME_FLOW*(OBX%INPUT_AREA(ABS(IOR))/OBX%UNDIVIDED_INPUT_AREA(ABS(IOR))) / OBX%FDS_AREA(ABS(IOR))
ENDIF
IF (OBST_INDEX>0 .AND. ABS(SF%MASS_FLUX_TOTAL)>=TWO_EPSILON_EB .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   OBX=>M%OBSTRUCTION(OBST_INDEX)
   WC%ONE_D%U_NORMAL_0 = SF%MASS_FLUX_TOTAL / RHOA * WC%ONE_D%AREA_ADJUST
ENDIF

! Do VENT-specific set-ups

T_ACTIVATE    = T_BEGIN-1._EB
WC%VENT_INDEX = VENT_INDEX

PROCESS_VENT: IF (WC%VENT_INDEX>0) THEN

   VT => M%VENTS(WC%VENT_INDEX)

   IF (.NOT.EVACUATION_ONLY(NM)) THEN
   WC%ONE_D%AREA_ADJUST = SF%AREA_MULTIPLIER * VT%INPUT_AREA/VT%FDS_AREA
   IF (WC%ONE_D%AREA_ADJUST<=TWO_EPSILON_EB) WC%ONE_D%AREA_ADJUST = 1._EB
   ENDIF

   ! Set the velocity at each surface cell

   WC%ONE_D%U_NORMAL_0 = SF%VEL

   IF (ABS(SF%VOLUME_FLOW)>TWO_EPSILON_EB)   THEN
      WC%ONE_D%U_NORMAL_0 = SF%VOLUME_FLOW*(VT%INPUT_AREA/VT%UNDIVIDED_INPUT_AREA)/VT%FDS_AREA
   ENDIF
   IF (ABS(SF%MASS_FLUX_TOTAL)>TWO_EPSILON_EB .AND. .NOT.EVACUATION_ONLY(NM)) THEN
      WC%ONE_D%U_NORMAL_0 = SF%MASS_FLUX_TOTAL/RHOA*WC%ONE_D%AREA_ADJUST
   ENDIF

   IF (SF%CONVERT_VOLUME_TO_MASS) THEN
      IF (ABS(WC%ONE_D%U_NORMAL_0)>TWO_EPSILON_EB) THEN
         ZZ_GET=0._EB
         ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,SF%MASS_FRACTION(1:N_TRACKED_SPECIES))
         CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F)
         SF%MASS_FLUX = -RHOA*(RSUM0/RSUM_F)*(TMPA/SF%TMP_FRONT)*SF%MASS_FRACTION*WC%ONE_D%U_NORMAL_0
         SF%SPECIES_BC_INDEX = SPECIFIED_MASS_FLUX
      ELSE
         CALL SHUTDOWN('ERROR: SURF: '//TRIM(SF%ID)//' must specify velocity boundary condition for conversion',&
                        PROCESS_0_ONLY=.FALSE.)
         IERR = 1
         RETURN
      ENDIF
   ENDIF

   ! Special velocity profiles

   PARABOLIC_IF: IF (SF%PROFILE==PARABOLIC_PROFILE) THEN
      SELECT CASE(ABS(IOR))
         CASE(1)
            IF (VT%RADIUS>0._EB) THEN
               RR = (M%YC(J)-VT%Y0)**2 + (M%ZC(K)-VT%Z0)**2
               WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*(VT%RADIUS**2-RR)/VT%RADIUS**2
            ELSE
               PY = 4._EB*(M%YC(J)-VT%Y1_ORIG)*(VT%Y2_ORIG-M%YC(J))/(VT%Y2_ORIG-VT%Y1_ORIG)**2
               PZ = 4._EB*(M%ZC(K)-VT%Z1_ORIG)*(VT%Z2_ORIG-M%ZC(K))/(VT%Z2_ORIG-VT%Z1_ORIG)**2
               WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*PY*PZ
            ENDIF
         CASE(2)
            IF (VT%RADIUS>0._EB) THEN
               RR = (M%XC(I)-VT%X0)**2 + (M%ZC(K)-VT%Z0)**2
               WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*(VT%RADIUS**2-RR)/VT%RADIUS**2
            ELSE
               PX = 4._EB*(M%XC(I)-VT%X1_ORIG)*(VT%X2_ORIG-M%XC(I))/(VT%X2_ORIG-VT%X1_ORIG)**2
               PZ = 4._EB*(M%ZC(K)-VT%Z1_ORIG)*(VT%Z2_ORIG-M%ZC(K))/(VT%Z2_ORIG-VT%Z1_ORIG)**2
               WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*PX*PZ
            ENDIF
         CASE(3)
            IF (VT%RADIUS>0._EB) THEN
               RR = (M%XC(I)-VT%X0)**2 + (M%YC(J)-VT%Y0)**2
               WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*(VT%RADIUS**2-RR)/VT%RADIUS**2
            ELSE
               PX = 4._EB*(M%XC(I)-VT%X1_ORIG)*(VT%X2_ORIG-M%XC(I))/(VT%X2_ORIG-VT%X1_ORIG)**2
               PY = 4._EB*(M%YC(J)-VT%Y1_ORIG)*(VT%Y2_ORIG-M%YC(J))/(VT%Y2_ORIG-VT%Y1_ORIG)**2
               IF (CYLINDRICAL) THEN
                  WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*PX
               ELSE
                  WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*PX*PY
               ENDIF
            ENDIF
      END SELECT
      IF (ABS(SF%VOLUME_FLOW)>=TWO_EPSILON_EB) THEN   ! Match desired volume flow
         IF (VT%RADIUS>0._EB) THEN
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*2._EB
         ELSE
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*9._EB/4._EB
         ENDIF
      ENDIF
   ENDIF PARABOLIC_IF

   IF (SF%PROFILE==BOUNDARY_LAYER_PROFILE) THEN

      ! Currently only set up for circular vents

      SELECT CASE(ABS(IOR))
         CASE(1)
            IF (VT%RADIUS>0._EB) THEN
               DELTA = VT%RADIUS - SQRT( VT%RADIUS**2*(2._EB*ABS(SF%VEL_BULK/SF%VEL)-1._EB) )
               R1 = VT%RADIUS - DELTA
               RR = SQRT( (M%YC(J)-VT%Y0)**2 + (M%ZC(K)-VT%Z0)**2 )
               IF (RR>R1 .AND. RR<=VT%RADIUS .AND. DELTA>TWO_EPSILON_EB) THEN
                  WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*(1._EB - ((RR-R1)/DELTA)**2 )
               ENDIF
            ENDIF
         CASE(2)
            IF (VT%RADIUS>0._EB) THEN
               DELTA = VT%RADIUS - SQRT( VT%RADIUS**2*(2._EB*ABS(SF%VEL_BULK/SF%VEL)-1._EB) )
               R1 = VT%RADIUS - DELTA
               RR = SQRT( (M%XC(I)-VT%X0)**2 + (M%ZC(K)-VT%Z0)**2 )
               IF (RR>R1 .AND. RR<=VT%RADIUS .AND. DELTA>TWO_EPSILON_EB) THEN
                  WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*(1._EB - ((RR-R1)/DELTA)**2 )
               ENDIF
            ENDIF
         CASE(3)
            IF (VT%RADIUS>0._EB) THEN
               DELTA = VT%RADIUS - SQRT( VT%RADIUS**2*(2._EB*ABS(SF%VEL_BULK/SF%VEL)-1._EB) )
               R1 = VT%RADIUS - DELTA
               RR = SQRT( (M%XC(I)-VT%X0)**2 + (M%YC(J)-VT%Y0)**2 )
               IF (RR>R1 .AND. RR<=VT%RADIUS .AND. DELTA>TWO_EPSILON_EB) THEN
                  WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*(1._EB - ((RR-R1)/DELTA)**2 )
               ENDIF
            ENDIF
      END SELECT
   ENDIF

   IF (SF%PROFILE==ATMOSPHERIC_PROFILE) THEN
      IF (M%ZC(K)<GROUND_LEVEL) THEN
         CALL SHUTDOWN('ERROR: SURF '//TRIM(SF%ID)//' cannot be applied below GROUND_LEVEL',PROCESS_0_ONLY=.FALSE.)
         IERR = 1
         RETURN
      ENDIF
      WC%ONE_D%U_NORMAL_0 =  WC%ONE_D%U_NORMAL_0*((M%ZC(K)-GROUND_LEVEL)/SF%Z0)**SF%PLE
   ENDIF

   IF (SF%PROFILE==RAMP_PROFILE) THEN
      SELECT CASE(ABS(IOR))
         CASE(1)
            IF (SF%RAMP_V_X/='null') THEN
               CALL SHUTDOWN('ERROR: RAMP_V_X assigned to SURF '//TRIM(SF%ID),PROCESS_0_ONLY=.FALSE.)
               IERR = 1
               RETURN
            ENDIF
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*EVALUATE_RAMP(M%YC(J),1._EB,SF%RAMP_INDEX(VELO_PROF_Y))
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*EVALUATE_RAMP(M%ZC(K),1._EB,SF%RAMP_INDEX(VELO_PROF_Z))
         CASE(2)
            IF (SF%RAMP_V_Y/='null') THEN
               CALL SHUTDOWN('ERROR: RAMP_V_Y assigned to SURF '//TRIM(SF%ID),PROCESS_0_ONLY=.FALSE.)
               IERR = 1
               RETURN
            ENDIF
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*EVALUATE_RAMP(M%XC(I),1._EB,SF%RAMP_INDEX(VELO_PROF_X))
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*EVALUATE_RAMP(M%ZC(K),1._EB,SF%RAMP_INDEX(VELO_PROF_Z))
         CASE(3)
            IF (SF%RAMP_V_Z/='null') THEN
               CALL SHUTDOWN('ERROR: RAMP_V_Z assigned to SURF '//TRIM(SF%ID),PROCESS_0_ONLY=.FALSE.)
               IERR = 1
               RETURN
            ENDIF
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*EVALUATE_RAMP(M%XC(I),1._EB,SF%RAMP_INDEX(VELO_PROF_X))
            WC%ONE_D%U_NORMAL_0 = WC%ONE_D%U_NORMAL_0*EVALUATE_RAMP(M%YC(J),1._EB,SF%RAMP_INDEX(VELO_PROF_Y))
      END SELECT
   ENDIF

   ! Check if fire spreads radially over this vent

   IF (VT%FIRE_SPREAD_RATE>0._EB) THEN
      DIST = SQRT((WC%X-VT%X0)**2 + (WC%Y-VT%Y0)**2 + (WC%Z-VT%Z0)**2)
      T_ACTIVATE = TT + DIST/VT%FIRE_SPREAD_RATE
   ENDIF

   ! Miscellaneous settings

   IF (.NOT.M%SOLID(ICG)) THEN
      IF (VT%BOUNDARY_TYPE==MIRROR_BOUNDARY) THEN
         WC%BOUNDARY_TYPE = MIRROR_BOUNDARY
         WC%SURF_INDEX    = MIRROR_SURF_INDEX
      ENDIF
      IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY) THEN
         WC%BOUNDARY_TYPE = OPEN_BOUNDARY
         WC%SURF_INDEX    = OPEN_SURF_INDEX
      ENDIF
   ENDIF

ENDIF PROCESS_VENT

! Check if fire spreads radially over this surface type

IF (SF%FIRE_SPREAD_RATE>0._EB) THEN
   DIST = SQRT((WC%X-SF%XYZ(1))**2 +(WC%Y-SF%XYZ(2))**2 +(WC%Z-SF%XYZ(3))**2)
   T_ACTIVATE = TT + DIST/SF%FIRE_SPREAD_RATE
ENDIF

! Set ignition time of each boundary cell

IF (T_ACTIVATE < T_BEGIN) THEN
   IF (SF%T_IGN==T_BEGIN) THEN
      WC%ONE_D%T_IGN = TT
   ELSE
      WC%ONE_D%T_IGN = SF%T_IGN
   ENDIF
ELSE
   WC%ONE_D%T_IGN = T_ACTIVATE
ENDIF

! Set correct initial value of temperature for RAMP_T

IF (ABS(WC%ONE_D%T_IGN-T_BEGIN) <= SPACING(WC%ONE_D%T_IGN) .AND. SF%RAMP_INDEX(TIME_TEMP)>=1) THEN
   TSI = TT
ELSE
   TSI = TT - WC%ONE_D%T_IGN
ENDIF

IF (SF%TMP_FRONT>0._EB) THEN
   WC%ONE_D%TMP_F = M%TMP_0(WC%ONE_D%KK) + &
                    EVALUATE_RAMP(TSI,SF%TAU(TIME_TEMP),SF%RAMP_INDEX(TIME_TEMP))*(SF%TMP_FRONT-M%TMP_0(WC%ONE_D%KK))
ELSE
   IF (SF%RAMP_T_I_INDEX < 0) WC%ONE_D%TMP_F = M%TMP_0(WC%ONE_D%KK)
ENDIF

! Record original boundary condition index for exterior wall cells that might get covered up

IF (OBST_INDEX==0) WC%SURF_INDEX_ORIG = SURF_INDEX_NEW

END SUBROUTINE INIT_WALL_CELL


SUBROUTINE OPEN_AND_CLOSE(T,NM)

! Check to see if a cell or OBSTruction is to be created or removed, or a VENT activated of deactivated

USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_STRINGS
USE CONTROL_VARIABLES, ONLY : CONTROL
USE DEVICE_VARIABLES, ONLY : DEVICE
USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
USE GEOMETRY_FUNCTIONS, ONLY: ASSIGN_HT3D_WALL_INDICES
REAL(EB), INTENT(IN) :: T
REAL(EB) :: T_TMP,TNOW
INTEGER  :: N,II,JJ,KK,IW,IC,VENT_INDEX,CVENT_INDEX
INTEGER, INTENT(IN) :: NM
LOGICAL :: CREATE_OBST,REMOVE_OBST,ACTIVATE_VENT,DEACTIVATE_VENT,ANY_REMOVE_OBST
CHARACTER(12) :: SV_LABEL
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (OBSTRUCTION_TYPE), POINTER :: OB

TNOW = CURRENT_TIME()

CALL POINT_TO_MESH(NM)

T_TMP = T
ANY_REMOVE_OBST = .FALSE.

! Check to see if an obstacle is to be removed or created

OBST_LOOP: DO N=1,N_OBST

   OB=>OBSTRUCTION(N)
   IF (.NOT. OB%REMOVABLE) CYCLE OBST_LOOP
   CREATE_OBST = .FALSE.
   REMOVE_OBST = .FALSE.

   ! Over-ride DEVICE/CONTROL logic

   CREATE_REMOVE_IF:IF (OB%CONSUMABLE .AND. OB%MASS<TWO_EPSILON_EB) THEN
      REMOVE_OBST = .TRUE.
   ELSE CREATE_REMOVE_IF
      SET_T_BEGIN_IF: IF (T<=T_BEGIN) THEN
         ! Set initial state of OBST
         HOLE_FILL_IF: IF (.NOT. OB%HOLE_FILLER) THEN
            !OBST is not a HOLE
            IF (OB%DEVC_INDEX > 0) THEN
               IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ELSEIF (OB%CTRL_INDEX > 0) THEN
               IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ENDIF
         ELSE HOLE_FILL_IF
            !OBST is a HOLE.  CREATE/REMOVE also depends on parent OBST.
            CHECK_PARENT: IF (OB%DEVC_INDEX_O > 0 .OR. OB%CTRL_INDEX_O > 0) THEN
               !Parent OBST controllable, check state and if parent OBST is hidden, do not fill hole.
               IF (OB%DEVC_INDEX_O > 0) THEN
                  IF (.NOT. DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ELSEIF(OB%CTRL_INDEX_O > 0) THEN
                  IF (.NOT. CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ENDIF
               !If parent OBST is visible, check to see if hole needs to be made.
               IF (.NOT. REMOVE_OBST) THEN
                  IF (OB%DEVC_INDEX > 0) THEN
                     IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                        REMOVE_OBST = .TRUE.
                     ELSE
                        CREATE_OBST = .TRUE.
                     ENDIF
                  ELSEIF (OB%CTRL_INDEX > 0) THEN
                     IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                        REMOVE_OBST = .TRUE.
                     ELSE
                        CREATE_OBST = .TRUE.
                     ENDIF
                  ENDIF
               ENDIF
            ELSE CHECK_PARENT
               !Parent OBST always present
               IF (OB%DEVC_INDEX > 0) THEN
                  IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ELSEIF (OB%CTRL_INDEX > 0) THEN
                  IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ENDIF
            ENDIF CHECK_PARENT

         ENDIF HOLE_FILL_IF
      ELSE SET_T_BEGIN_IF
      ! Decide if a DEVICE/CONTROL action is needed
         HOLE_FILL_IF_2: IF (.NOT. OB%HOLE_FILLER) THEN
            !OBST is not a HOLE
            IF (OB%DEVC_INDEX > 0) THEN
               IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
               IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ELSEIF (OB%CTRL_INDEX > 0) THEN
               IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
               IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ENDIF
         ELSE HOLE_FILL_IF_2
            !OBST is a HOLE.  CREATE/REMOVE also depends on parent OBST.
            CHECK_PARENT_2: IF (OB%DEVC_INDEX_O > 0 .OR. OB%CTRL_INDEX_O > 0) THEN
               !Parent OBST controllable, check state and if parent OBST is hidden, do not fill hole.
               IF (OB%DEVC_INDEX_O > 0) THEN
                  IF (.NOT. DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ELSEIF(OB%CTRL_INDEX_O > 0) THEN
                  IF (.NOT. CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ENDIF
               !If parent OBST is visible, check to see if hole needs to be made.
               IF (.NOT. REMOVE_OBST) THEN
                  IF (OB%DEVC_INDEX > 0) THEN
                     IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) THEN
                        IF (OB%DEVC_INDEX_O > 0 .AND. .NOT. DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                           IF (DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE .NEQV. DEVICE(OB%DEVC_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ELSEIF(OB%CTRL_INDEX_O > 0 .AND. .NOT. DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                           IF (CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE .NEQV. CONTROL(OB%CTRL_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ENDIF
                     ELSE
                        IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                           REMOVE_OBST = .TRUE.
                        ELSE
                           CREATE_OBST = .TRUE.
                        ENDIF
                     ENDIF
                  ELSEIF (OB%CTRL_INDEX > 0) THEN
                     IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) THEN
                        IF (OB%DEVC_INDEX_O > 0 .AND. .NOT. CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                           IF (DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE .NEQV. DEVICE(OB%DEVC_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ELSEIF(OB%CTRL_INDEX_O > 0 .AND. .NOT. CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                           IF (CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE .NEQV. CONTROL(OB%CTRL_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ENDIF
                     ELSE
                        IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                           REMOVE_OBST = .TRUE.
                        ELSE
                           CREATE_OBST = .TRUE.
                        ENDIF
                     ENDIF
                  ENDIF
               ENDIF
            ELSE CHECK_PARENT_2
               !Parent OBST not controllable and is always present
               IF (OB%DEVC_INDEX > 0) THEN
                  IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
                  IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ELSEIF (OB%CTRL_INDEX > 0) THEN
                  IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
                  IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ENDIF
            ENDIF CHECK_PARENT_2

         ENDIF HOLE_FILL_IF_2

      ENDIF SET_T_BEGIN_IF
   ENDIF CREATE_REMOVE_IF

   SV_LABEL  = 'null'

   IF (CREATE_OBST .AND. OB%HIDDEN) THEN
      OB%HIDDEN = .FALSE.
      SV_LABEL  = 'SHOW_OBST'
      CALL CREATE_OR_REMOVE_OBST(NM,OB%I1,OB%I2,OB%J1,OB%J2,OB%K1,OB%K2,1,N)
   ENDIF

   IF (REMOVE_OBST .AND. (.NOT. OB%HIDDEN)) THEN
      OB%HIDDEN = .TRUE.
      SV_LABEL  = 'HIDE_OBST'
      CALL CREATE_OR_REMOVE_OBST(NM,OB%I1,OB%I2,OB%J1,OB%J2,OB%K1,OB%K2,0,N)
   ENDIF

   ! The evacuation flow field calculation is done before T_BEGIN

   IF (OB%DEVC_INDEX>0) THEN
      IF (EVACUATION_SKIP(NM) .AND. DEVICE(OB%DEVC_INDEX)%QUANTITY(1)=='TIME' .AND. DEVICE(OB%DEVC_INDEX)%SETPOINT<=T_BEGIN) &
         T_TMP = T - EVAC_DT_FLOWFIELD*EVAC_TIME_ITERATIONS
   ENDIF

   ! Write a message to the Smokeview .smv file that the obstruction has been created or removed

   IF (SV_LABEL /= 'null') THEN
      IF (N_STRINGS+2>N_STRINGS_MAX) THEN
         CALL RE_ALLOCATE_STRINGS(NM)
         STRING => MESHES(NM)%STRING
      ENDIF
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(I6,F14.6)') N,T_BEGIN+(T_TMP-T_BEGIN)*TIME_SHRINK_FACTOR
   ENDIF

   ! If any REMOVE_OBST store

   IF (REMOVE_OBST) ANY_REMOVE_OBST=.TRUE.

ENDDO OBST_LOOP

! Reassign wall indices in PYRO3D

IF (ANY_REMOVE_OBST .AND. SOLID_PYRO3D .AND. .NOT.SOLID_MT3D) CALL ASSIGN_HT3D_WALL_INDICES(NM)

! Check to see if a vent should be activated or deactivated

VENT_INDEX  = 0
CVENT_INDEX = 0

VENT_LOOP: DO N=1,N_VENT
   VT => VENTS(N)

   IF (VT%RADIUS>0._EB) THEN
      CVENT_INDEX = CVENT_INDEX + 1
   ELSE
      VENT_INDEX = VENT_INDEX + 1
   ENDIF

   ACTIVATE_VENT   = .FALSE.
   DEACTIVATE_VENT = .FALSE.

   ! Over-ride DEVICE/CONTROL logic

   IF (.NOT.VT%ACTIVATED .AND. T<=T_BEGIN) DEACTIVATE_VENT = .TRUE.

   ! If the VENT is tied to a specific OBST, and the OBST is HIDDEN (not HIDDEN), and the VENT is activated (not activated),
   ! deactivate (activate) the vent.

   IF (VT%OBST_INDEX>0 .AND. OBSTRUCTION(VT%OBST_INDEX)%HIDDEN .AND. VT%ACTIVATED) THEN
      VT%ACTIVATED = .FALSE.
      DEACTIVATE_VENT = .TRUE.
   ENDIF

   IF (VT%OBST_INDEX>0 .AND. .NOT.OBSTRUCTION(VT%OBST_INDEX)%HIDDEN .AND. .NOT.VT%ACTIVATED) THEN
      VT%ACTIVATED = .TRUE.
      ACTIVATE_VENT = .TRUE.
   ENDIF

   ! Decide if a VENT is to activate or de-activate based on a DEVICE or CONTROLLER

   IF (.NOT.ACTIVATE_VENT .AND. .NOT.DEACTIVATE_VENT) THEN
      IF (VT%DEVC_INDEX > 0) THEN
         IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(VT%DEVC_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
         IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE) THEN
            ACTIVATE_VENT   = .TRUE.
         ELSE
            DEACTIVATE_VENT = .TRUE.
         ENDIF
      ELSEIF (VT%CTRL_INDEX > 0) THEN
         IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(VT%CTRL_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
         IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE) THEN
            ACTIVATE_VENT   = .TRUE.
         ELSE
            DEACTIVATE_VENT = .TRUE.
         ENDIF
      ENDIF
   ENDIF

   IF (.NOT.ACTIVATE_VENT .AND. .NOT.DEACTIVATE_VENT) CYCLE VENT_LOOP

   ! Find the wall indices (IW) for the vent and set the activation time (ONE_D%T_IGN) for each one

   DO KK=VT%K1+1,MAX(VT%K1+1,VT%K2)
      DO JJ=VT%J1+1,MAX(VT%J1+1,VT%J2)
         DO II=VT%I1+1,MAX(VT%I1+1,VT%I2)
            SELECT CASE(VT%IOR)
               CASE(1:)
                  IC = CELL_INDEX(II,JJ,KK)
               CASE(-1)
                  IC = CELL_INDEX(II-1,JJ,KK)
               CASE(-2)
                  IC = CELL_INDEX(II,JJ-1,KK)
               CASE(-3)
                  IC = CELL_INDEX(II,JJ,KK-1)
            END SELECT
            IW = WALL_INDEX(IC,-VT%IOR)
            IF (IW==0) CYCLE

            IF (ACTIVATE_VENT) THEN
               IF (VT%FIRE_SPREAD_RATE>0._EB) THEN
                  WALL(IW)%ONE_D%T_IGN = T + SQRT((WALL(IW)%X-VT%X0)**2 +(WALL(IW)%Y-VT%Y0)**2 +&
                                              (WALL(IW)%Z-VT%Z0)**2)/VT%FIRE_SPREAD_RATE
               ELSE
                  WALL(IW)%ONE_D%T_IGN = T
               ENDIF
            ELSE
               WALL(IW)%ONE_D%T_IGN = 1.E6_EB
            ENDIF
         ENDDO
      ENDDO
   ENDDO

   ! Write message to .smv file

   IF (VT%RADIUS<0._EB) THEN
      IF (ACTIVATE_VENT)   SV_LABEL = 'OPEN_VENT'
      IF (DEACTIVATE_VENT) SV_LABEL = 'CLOSE_VENT'
   ELSE
      IF (ACTIVATE_VENT)   SV_LABEL = 'OPEN_CVENT'
      IF (DEACTIVATE_VENT) SV_LABEL = 'CLOSE_CVENT'
   ENDIF

   IF (N_STRINGS+2>N_STRINGS_MAX) THEN
      CALL RE_ALLOCATE_STRINGS(NM)
      STRING => MESHES(NM)%STRING
   ENDIF
   N_STRINGS = N_STRINGS + 1
   WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
   N_STRINGS = N_STRINGS + 1
   IF (VT%RADIUS>0._EB) WRITE(STRING(N_STRINGS),'(I6,F10.2)') CVENT_INDEX,T
   IF (VT%RADIUS<0._EB) WRITE(STRING(N_STRINGS),'(I6,F10.2)')  VENT_INDEX,T

ENDDO VENT_LOOP

T_USED(6) = T_USED(6) + CURRENT_TIME() - TNOW
END SUBROUTINE OPEN_AND_CLOSE


SUBROUTINE CREATE_OR_REMOVE_OBST(NM,I1,I2,J1,J2,K1,K2,CR_INDEX,OBST_INDEX)

! Create or remove the obstruction whose NODES (not cells) are given by I1, I2, etc.

USE GEOMETRY_FUNCTIONS, ONLY : BLOCK_CELL
INTEGER :: I1,I2,J1,J2,K1,K2,I,J,K
INTEGER, INTENT(IN) :: NM,CR_INDEX,OBST_INDEX
LOGICAL :: CREATE,REMOVE

CALL POINT_TO_MESH(NM)

! Indicate whether to create or remove the obstruction.

OBST_CREATED_OR_REMOVED = .TRUE.
REMOVE = .FALSE.
CREATE = .FALSE.
IF (CR_INDEX==0) REMOVE = .TRUE.
IF (CR_INDEX==1) CREATE = .TRUE.
IF (REMOVE) OBSTRUCTION(OBST_INDEX)%SCHEDULED_FOR_REMOVAL  = .TRUE.
IF (CREATE) OBSTRUCTION(OBST_INDEX)%SCHEDULED_FOR_CREATION = .TRUE.

! Blank or unblank cells that make up the OBSTruction

IF (I1/=I2 .AND. J1/=J2 .AND. K1/=K2) CALL BLOCK_CELL(NM,I1+1,I2,J1+1,J2,K1+1,K2,CR_INDEX,OBST_INDEX)

! If the OBSTruction is to be removed, set density and mass fractions to ambient value

IF (REMOVE) THEN
   DO K=K1+1,K2
      DO J=J1+1,J2
         DO I=I1+1,I2
            RHOS(I,J,K) = RHOA
            RHO(I,J,K)  = RHOA
            IF (SOLID_PHASE_ONLY) TMP(I,J,K)  = TMPA
            ZZ(I,J,K,1:N_TRACKED_SPECIES)  = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
            ZZS(I,J,K,1:N_TRACKED_SPECIES) = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
         ENDDO
      ENDDO
   ENDDO
ENDIF

END SUBROUTINE CREATE_OR_REMOVE_OBST


!> \brief Re-assign wall boundaries after create or removal of obstructions

SUBROUTINE REASSIGN_WALL_CELLS(T,NM)

USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: T
INTEGER :: N,I1,I2,J1,J2,K1,K2,I,J,K,IW,ICG,IC,OBST_INDEX,NOM,IIO,JJO,KKO
REAL(EB) :: TNOW
LOGICAL :: CREATE,REMOVE
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (WALL_TYPE), POINTER :: WC
TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE (MESH_TYPE), POINTER :: MM

TNOW = CURRENT_TIME()

CALL POINT_TO_MESH(NM)

DO IW=1,N_EXTERNAL_WALL_CELLS
   WC => WALL(IW)
   EWC => EXTERNAL_WALL(IW)
   NOM = EWC%NOM
   IF (NOM==0) CYCLE
   MM  => MESHES(NOM)
   IIO = EWC%IIO_MIN
   JJO = EWC%JJO_MIN
   KKO = EWC%KKO_MIN
   IF (WC%OBST_INDEX==0 .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. .NOT.MM%SOLID(MM%CELL_INDEX(IIO,JJO,KKO))) THEN
      IC  = CELL_INDEX(WC%ONE_D%II,WC%ONE_D%JJ,WC%ONE_D%KK)
      ICG = CELL_INDEX(WC%ONE_D%IIG,WC%ONE_D%JJG,WC%ONE_D%KKG)
      IF (SOLID(ICG)) CYCLE
      OBST_INDEX=0
      CALL GET_BOUNDARY_TYPE
   ENDIF
ENDDO

! Loop over all obstructions in the current mesh and initialize newly exposed or covered wall cell faces

OBSTRUCTION_LOOP: DO N=1,N_OBST

OB => OBSTRUCTION(N)
OBST_INDEX=N

IF (.NOT.OB%SCHEDULED_FOR_REMOVAL .AND. .NOT.OB%SCHEDULED_FOR_CREATION) CYCLE OBSTRUCTION_LOOP

REMOVE = .FALSE. ; CREATE = .FALSE.
IF (OB%SCHEDULED_FOR_REMOVAL) THEN
   REMOVE = .TRUE.
   OB%SCHEDULED_FOR_REMOVAL = .FALSE.
ENDIF
IF (OB%SCHEDULED_FOR_CREATION) THEN
   CREATE = .TRUE.
   OB%SCHEDULED_FOR_CREATION = .FALSE.
ENDIF

I1 = OB%I1 ; I2 = OB%I2 ; J1 = OB%J1 ; J2 = OB%J2 ; K1 = OB%K1 ; K2 = OB%K2

DO K=K1+1,K2
   DO J=J1+1,J2
      IC  = CELL_INDEX(I1+1,J,K)
      ICG = CELL_INDEX(I1  ,J,K)
      IW  = WALL_INDEX(ICG,1)
      IF (I1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I1  ,J,K)
      ICG = CELL_INDEX(I1+1,J,K)
      IW  = WALL_INDEX(ICG,-1)
      IF (I1<IBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I2  ,J,K)
      ICG = CELL_INDEX(I2+1,J,K)
      IW  = WALL_INDEX(ICG,-1)
      IF (I2<IBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I2+1,J,K)
      ICG = CELL_INDEX(I2  ,J,K)
      IW  = WALL_INDEX(ICG, 1)
      IF (I2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Process the y boundaries of the OBSTruction

DO K=K1+1,K2
   DO I=I1+1,I2
      IC  = CELL_INDEX(I,J1+1,K)
      ICG = CELL_INDEX(I,J1  ,K)
      IW  = WALL_INDEX(ICG, 2)
      IF (J1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J1  ,K)
      ICG = CELL_INDEX(I,J1+1,K)
      IW  = WALL_INDEX(ICG,-2)
      IF (J1<JBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J2  ,K)
      ICG = CELL_INDEX(I,J2+1,K)
      IW  = WALL_INDEX(ICG,-2)
      IF (J2<JBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J2+1,K)
      ICG = CELL_INDEX(I,J2  ,K)
      IW  = WALL_INDEX(ICG, 2)
      IF (J2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Process the z boundaries of the OBSTruction

DO J=J1+1,J2
   DO I=I1+1,I2
      IC  = CELL_INDEX(I,J,K1+1)
      ICG = CELL_INDEX(I,J,K1  )
      IW  = WALL_INDEX(ICG, 3)
      IF (K1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K1  )
      ICG = CELL_INDEX(I,J,K1+1)
      IW  = WALL_INDEX(ICG,-3)
      IF (K1<KBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K2  )
      ICG = CELL_INDEX(I,J,K2+1)
      IW  = WALL_INDEX(ICG,-3)
      IF (K2<KBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K2+1)
      ICG = CELL_INDEX(I,J,K2  )
      IW  = WALL_INDEX(ICG, 3)
      IF (K2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Nullify block edges on blockage that is to be removed

DO K=K1,K2
   DO J=J1,J2
      IF (J>J1) CALL REDEFINE_EDGE(I1,J,K,2)
      IF (J>J1) CALL REDEFINE_EDGE(I2,J,K,2)
      IF (K>K1) CALL REDEFINE_EDGE(I1,J,K,3)
      IF (K>K1) CALL REDEFINE_EDGE(I2,J,K,3)
   ENDDO
ENDDO

DO K=K1,K2
   DO I=I1,I2
      IF (I>I1) CALL REDEFINE_EDGE(I,J1,K,1)
      IF (I>I1) CALL REDEFINE_EDGE(I,J2,K,1)
      IF (K>K1) CALL REDEFINE_EDGE(I,J1,K,3)
      IF (K>K1) CALL REDEFINE_EDGE(I,J2,K,3)
   ENDDO
ENDDO

DO J=J1,J2
   DO I=I1,I2
      IF (I>I1) CALL REDEFINE_EDGE(I,J,K1,1)
      IF (I>I1) CALL REDEFINE_EDGE(I,J,K2,1)
      IF (J>J1) CALL REDEFINE_EDGE(I,J,K1,2)
      IF (J>J1) CALL REDEFINE_EDGE(I,J,K2,2)
   ENDDO
ENDDO

ENDDO OBSTRUCTION_LOOP

T_USED(6) = T_USED(6) + CURRENT_TIME() - TNOW
CONTAINS

SUBROUTINE GET_BOUNDARY_TYPE

INTEGER :: IOR,IIG,JJG,KKG,IW_OLD,IERR,PRESSURE_BC_INDEX

IF (IW==0) RETURN

WC => WALL(IW)

IF (IW<=N_EXTERNAL_WALL_CELLS .AND. REMOVE) THEN
   WC%BOUNDARY_TYPE = SOLID_BOUNDARY
   SOLID(IC)        = .TRUE.
   IF (WC%SURF_INDEX_ORIG==MIRROR_SURF_INDEX) THEN
      WC%BOUNDARY_TYPE = MIRROR_BOUNDARY
      WC%SURF_INDEX    = MIRROR_SURF_INDEX
      SOLID(IC)        = .TRUE.
      RETURN
   ENDIF
   IF (WC%SURF_INDEX_ORIG==OPEN_SURF_INDEX) THEN
      WC%BOUNDARY_TYPE = OPEN_BOUNDARY
      WC%SURF_INDEX    = OPEN_SURF_INDEX
      SOLID(IC)        = .FALSE.
      RETURN
   ENDIF
   IF (WC%SURF_INDEX_ORIG==INTERPOLATED_SURF_INDEX) THEN
      WC%BOUNDARY_TYPE = INTERPOLATED_BOUNDARY
      WC%SURF_INDEX    = INTERPOLATED_SURF_INDEX
      SOLID(IC)        = .FALSE.
      RETURN
   ENDIF
ENDIF

IF (IW>N_EXTERNAL_WALL_CELLS) THEN
   IF (WC%OBST_INDEX>0 .AND.      OBSTRUCTION(WC%OBST_INDEX)%HIDDEN .AND. .NOT.SOLID(IC) ) WC%BOUNDARY_TYPE = NULL_BOUNDARY
   IF (WC%OBST_INDEX>0 .AND. .NOT.OBSTRUCTION(WC%OBST_INDEX)%HIDDEN .AND. .NOT.SOLID(ICG)) WC%BOUNDARY_TYPE = SOLID_BOUNDARY
   IF (SOLID(ICG)) WC%BOUNDARY_TYPE = NULL_BOUNDARY
ENDIF

IF (CREATE) THEN
   IF (SOLID(ICG)) THEN
      WC%BOUNDARY_TYPE = NULL_BOUNDARY
   ELSE
      WC%BOUNDARY_TYPE = SOLID_BOUNDARY
      IF (WC%ONE_D%T_IGN<T) WC%ONE_D%T_IGN = T
   ENDIF
ENDIF

! Initialize solid properties included in the WALL derived type variable

IF (WC%BOUNDARY_TYPE/=NULL_BOUNDARY) THEN
   IF (IW<=N_EXTERNAL_WALL_CELLS) THEN
      WC%OBST_INDEX = 0
      WC%SURF_INDEX = WC%SURF_INDEX_ORIG
   ENDIF
   IF (OBST_INDEX_C(IC)>0) THEN
      WC%OBST_INDEX = OBST_INDEX_C(IC)
      WC%SURF_INDEX = OBSTRUCTION(WC%OBST_INDEX)%SURF_INDEX(WC%ONE_D%IOR)
   ELSEIF (CREATE .AND. OBST_INDEX>0) THEN
      WC%OBST_INDEX = OBST_INDEX
      WC%SURF_INDEX = OBSTRUCTION(WC%OBST_INDEX)%SURF_INDEX(WC%ONE_D%IOR)
   ENDIF
   PRESSURE_BC_INDEX = WC%PRESSURE_BC_INDEX  ! Save this parameter and restore it after the call to INIT_WALL_CELL
   CALL INIT_WALL_CELL(NM,WC%ONE_D%II,WC%ONE_D%JJ,WC%ONE_D%KK,WC%OBST_INDEX,IW,WC%ONE_D%IOR,WC%SURF_INDEX,IERR,T)
   WC => WALL(IW)
   WC%PRESSURE_BC_INDEX = PRESSURE_BC_INDEX
ENDIF

! Special case: A BURNed_AWAY obstruction exposes a surface that also burns, in which case the surface is to ignite immediately.

IF (REMOVE .AND. SURFACE(WC%SURF_INDEX)%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED) THEN
   IIG = WC%ONE_D%IIG
   JJG = WC%ONE_D%JJG
   KKG = WC%ONE_D%KKG
   IOR = WC%ONE_D%IOR
   IW_OLD = 0
   IC     = 0
   SELECT CASE(IOR)
      CASE(-1)
         IF (IIG>1)    IC = CELL_INDEX(IIG-1,JJG,KKG)
      CASE( 1)
         IF (IIG<IBAR) IC = CELL_INDEX(IIG+1,JJG,KKG)
      CASE(-2)
         IF (JJG>1)    IC = CELL_INDEX(IIG,JJG-1,KKG)
      CASE( 2)
         IF (JJG<JBAR) IC = CELL_INDEX(IIG,JJG+1,KKG)
      CASE(-3)
         IF (KKG>1)    IC = CELL_INDEX(IIG,JJG,KKG-1)
      CASE( 3)
         IF (KKG<KBAR) IC = CELL_INDEX(IIG,JJG,KKG+1)
   END SELECT
   IW_OLD = WALL_INDEX(IC,-IOR)
   IF (IW_OLD>0) THEN
      IF (WALL(IW_OLD)%SURF_INDEX==WC%SURF_INDEX) WC%ONE_D%T_IGN =  WALL(IW_OLD)%ONE_D%T_IGN
   ENDIF
ENDIF

END SUBROUTINE GET_BOUNDARY_TYPE


SUBROUTINE REDEFINE_EDGE(II,JJ,KK,IEC)

! Change a few properties of the EDGEs that have been exposed or covered up by a blockage

INTEGER :: IE,II,JJ,KK,IEC

SELECT CASE(IEC)
   CASE(1)
      IE = EDGE_INDEX( 4,CELL_INDEX(II,JJ,KK))
   CASE(2)
      IE = EDGE_INDEX( 8,CELL_INDEX(II,JJ,KK))
   CASE(3)
      IE = EDGE_INDEX(12,CELL_INDEX(II,JJ,KK))
END SELECT

END SUBROUTINE REDEFINE_EDGE

END SUBROUTINE REASSIGN_WALL_CELLS


SUBROUTINE DEFINE_EDGE(II,JJ,KK,IOR,IEC,NM,OBST_INDEX,IERR)

! Set up edge arrays for velocity boundary conditions

INTEGER, INTENT(IN) :: II,JJ,KK,IOR,IEC,NM
INTEGER :: NOM,ICMM,ICMP,ICPM,ICPP,OBST_INDEX,IE,IW,IIO,JJO,KKO,IW1,IW2,IERR
REAL(EB) :: XI,YJ,ZK
TYPE (MESH_TYPE), POINTER :: M,MM
TYPE (OBSTRUCTION_TYPE), POINTER :: OB

M => MESHES(NM)
IF (OBST_INDEX>0) OB=>OBSTRUCTION(OBST_INDEX)

! Find the wall cells on each side of the edge

IW1 = -1
IW2 = -1

EDGE_DIRECTION_1: SELECT CASE(IEC)
   CASE(1) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-2)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ,KK)  ,2)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ,KK+1),2)
         CASE( 2)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK)  ,-2)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK+1),-2)
         CASE(-3)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ  ,KK),3)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK),3)
         CASE( 3)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ  ,KK+1),-3)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK+1),-3)
      END SELECT
   CASE(2) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-1)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ,KK)  ,1)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ,KK+1),1)
         CASE( 1)
            IW1 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK)  ,-1)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK+1),-1)
         CASE(-3)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ,KK),3)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK),3)
         CASE( 3)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ,KK+1),-3)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK+1),-3)
      END SELECT
   CASE(3) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-1)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ  ,KK),1)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK),1)
         CASE( 1)
            IW1 = WALL_INDEX(CELL_INDEX(II+1,JJ  ,KK),-1)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ+1,KK),-1)
         CASE(-2)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ,KK),2)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK),2)
         CASE( 2)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ+1,KK),-2)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ+1,KK),-2)
      END SELECT
END SELECT EDGE_DIRECTION_1

! Decide what to do based on whether or not adjacent tiles exist

IF (IW1==-1 .OR. IW2==-1) THEN
   WRITE(LU_ERR,'(A,I2,A,3I3)') 'ERROR: Edge initialization failed; Mesh: ',NM,', Cell: ',II,JJ,KK
   STOP_STATUS = SETUP_STOP
   IERR = 1
   RETURN
ENDIF

IF (IW1==0 .AND. IW2==0) RETURN
IF (IW1> 0 .AND. IW2==0) IW = IW1
IF (IW1==0 .AND. IW2> 0) IW = IW2
IF (IW1> 0 .AND. IW2> 0) THEN
   IW = IW2
   IF (IW1<=N_EXTERNAL_WALL_CELLS) THEN
      IF (EXTERNAL_WALL(IW1)%NOM>0) IW = IW1
   ENDIF
   IF (IW2<=N_EXTERNAL_WALL_CELLS) THEN
      IF (EXTERNAL_WALL(IW2)%NOM>0) IW = IW2
   ENDIF
ENDIF

! Assign the Index of the Edge (IE) and add to the list

ICMM = CELL_INDEX(II,JJ,KK)
SELECT CASE(IEC)
   CASE(1)
      IE = EDGE_INDEX( 4,ICMM)
   CASE(2)
      IE = EDGE_INDEX( 8,ICMM)
   CASE(3)
      IE = EDGE_INDEX(12,ICMM)
END SELECT

IF (IE==0) THEN
   N_EDGES = N_EDGES + 1
   IE = N_EDGES
ENDIF

! Determine the wall index of the adjacent wall tile

NOM = 0
IIO = 0
JJO = 0
KKO = 0

IF (IW<=N_EXTERNAL_WALL_CELLS) THEN
   IF (EXTERNAL_WALL(IW)%NOM>0) THEN
      NOM = EXTERNAL_WALL(IW)%NOM
      IIO = EXTERNAL_WALL(IW)%IIO_MIN
      JJO = EXTERNAL_WALL(IW)%JJO_MIN
      KKO = EXTERNAL_WALL(IW)%KKO_MIN
   ENDIF
ENDIF

! Fill up array IJKE with edge parameters

IJKE(1,IE) = II
IJKE(2,IE) = JJ
IJKE(3,IE) = KK
IJKE(4,IE) = IEC

! Fill in EDGE_INDEX and the rest of IJKE

EDGE_DIRECTION_2: SELECT CASE(IEC)

   CASE (1) EDGE_DIRECTION_2

      ICPM = CELL_INDEX(II,JJ+1,KK)
      ICPP = CELL_INDEX(II,JJ+1,KK+1)
      ICMP = CELL_INDEX(II,JJ,KK+1)
      IJKE(5,IE) = ICMM
      IJKE(6,IE) = ICPM
      IJKE(7,IE) = ICMP
      IJKE(8,IE) = ICPP
      EDGE_INDEX(1,ICPP) = IE
      EDGE_INDEX(2,ICMP) = IE
      EDGE_INDEX(3,ICPM) = IE
      EDGE_INDEX(4,ICMM) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE(2)
               IF (IOR>0) IJKE( 9,IE) = -NOM
               IF (IOR<0) IJKE( 9,IE) =  NOM
               IJKE(10,IE) = IIO
               IJKE(11,IE) = JJO
               MM => MESHES(NOM)
               ZK  = MIN( REAL(MM%KBAR,EB)+ALMOST_ONE , MM%CELLSK(NINT((Z(KK)-MM%ZS)*MM%RDZINT))+1._EB )
               KKO = MAX(1,FLOOR(ZK))
               M%EDGE_INTERPOLATION_FACTOR(IE,1) = ZK-KKO
               IJKE(12,IE) = KKO
            CASE(3)
               IF (IOR>0) IJKE(13,IE) = -NOM
               IF (IOR<0) IJKE(13,IE) =  NOM
               IJKE(14,IE) = IIO
               MM => MESHES(NOM)
               YJ  = MIN( REAL(MM%JBAR,EB)+ALMOST_ONE , MM%CELLSJ(NINT((Y(JJ)-MM%YS)*MM%RDYINT))+1._EB )
               JJO = MAX(1,FLOOR(YJ))
               M%EDGE_INTERPOLATION_FACTOR(IE,2) = YJ-JJO
               IJKE(15,IE) = JJO
               IJKE(16,IE) = KKO
         END SELECT
      ENDIF

   CASE (2) EDGE_DIRECTION_2

      ICMP = CELL_INDEX(II+1,JJ,KK)
      ICPP = CELL_INDEX(II+1,JJ,KK+1)
      ICPM = CELL_INDEX(II,JJ,KK+1)
      IJKE(5,IE) = ICMM
      IJKE(6,IE) = ICPM
      IJKE(7,IE) = ICMP
      IJKE(8,IE) = ICPP
      EDGE_INDEX(5,ICPP) = IE
      EDGE_INDEX(6,ICPM) = IE
      EDGE_INDEX(7,ICMP) = IE
      EDGE_INDEX(8,ICMM) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE( 1)
               IF (IOR>0) IJKE(13,IE) = -NOM
               IF (IOR<0) IJKE(13,IE) =  NOM
               IJKE(14,IE) = IIO
               IJKE(15,IE) = JJO
               MM => MESHES(NOM)
               ZK  = MIN( REAL(MM%KBAR,EB)+ALMOST_ONE , MM%CELLSK(NINT((Z(KK)-MM%ZS)*MM%RDZINT))+1._EB )
               KKO = MAX(1,FLOOR(ZK))
               M%EDGE_INTERPOLATION_FACTOR(IE,2) = ZK-KKO
               IJKE(16,IE) = KKO
            CASE( 3)
               IF (IOR>0) IJKE( 9,IE) = -NOM
               IF (IOR<0) IJKE( 9,IE) =  NOM
               MM => MESHES(NOM)
               XI  = MIN( REAL(MM%IBAR,EB)+ALMOST_ONE , MM%CELLSI(NINT((X(II)-MM%XS)*MM%RDXINT))+1._EB )
               IIO = MAX(1,FLOOR(XI))
               M%EDGE_INTERPOLATION_FACTOR(IE,1) = XI-IIO
               IJKE(10,IE) = IIO
               IJKE(11,IE) = JJO
               IJKE(12,IE) = KKO
         END SELECT
      ENDIF

   CASE (3) EDGE_DIRECTION_2

      ICPM = CELL_INDEX(II+1,JJ,KK)
      ICPP = CELL_INDEX(II+1,JJ+1,KK)
      ICMP = CELL_INDEX(II,JJ+1,KK)
      IJKE(5,IE) = ICMM
      IJKE(6,IE) = ICPM
      IJKE(7,IE) = ICMP
      IJKE(8,IE) = ICPP
      EDGE_INDEX( 9,ICPP) = IE
      EDGE_INDEX(10,ICMP) = IE
      EDGE_INDEX(11,ICPM) = IE
      EDGE_INDEX(12,ICMM) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE( 1)
               IF (IOR>0) IJKE( 9,IE) = -NOM
               IF (IOR<0) IJKE( 9,IE) =  NOM
               IJKE(10,IE) = IIO
               MM => MESHES(NOM)
               YJ  = MIN( REAL(MM%JBAR,EB)+ALMOST_ONE , MM%CELLSJ(NINT((Y(JJ)-MM%YS)*MM%RDYINT))+1._EB )
               JJO = MAX(1,FLOOR(YJ))
               M%EDGE_INTERPOLATION_FACTOR(IE,1) = YJ-JJO
               IJKE(11,IE) = JJO
               IJKE(12,IE) = KKO
            CASE( 2)
               IF (IOR>0) IJKE(13,IE) = -NOM
               IF (IOR<0) IJKE(13,IE) =  NOM
               MM => MESHES(NOM)
               XI  = MIN( REAL(MM%IBAR,EB)+ALMOST_ONE , MM%CELLSI(NINT((X(II)-MM%XS)*MM%RDXINT))+1._EB )
               IIO = MAX(1,FLOOR(XI))
               M%EDGE_INTERPOLATION_FACTOR(IE,2) = XI-IIO
               IJKE(14,IE) = IIO
               IJKE(15,IE) = JJO
               IJKE(16,IE) = KKO
         END SELECT
      ENDIF

END SELECT EDGE_DIRECTION_2

END SUBROUTINE DEFINE_EDGE


SUBROUTINE INITIAL_NOISE(NM)

! Generate random noise at the start of the simulation

REAL     :: RN2
REAL(EB) :: RN
INTEGER  :: I,J,K,SIZE_RND,IZERO
INTEGER, DIMENSION(:), ALLOCATABLE :: SEED_RND
INTEGER, INTENT(IN) :: NM

IF (EVACUATION_ONLY(NM)) RETURN

! Waste a few calls to RANDOM_NUMBER to avoid generating the exact same sequence on each mesh

CALL RANDOM_SEED(SIZE=SIZE_RND)
ALLOCATE(SEED_RND(SIZE_RND),STAT=IZERO)
CALL CHKMEMERR('INITIAL_NOISE','SEED_RND',IZERO)
SEED_RND = 2819 * 13*NM
CALL RANDOM_SEED(PUT=SEED_RND)
DEALLOCATE(SEED_RND)

DO I=1,NM
   CALL RANDOM_NUMBER(RN2)
ENDDO

! Point to local mesh variables

CALL POINT_TO_MESH(NM)

! Add random vorticity to cells that are not bounding solid surfaces

DO K=1,KBM1
   DO J=1,JBM1
      DO I=1,IBAR
      !  IF (WALL_INDEX(CELL_INDEX(I,J  ,K), 3)>0 .OR. WALL_INDEX(CELL_INDEX(I,J+1,K), 3)>0 .OR. &
      !      WALL_INDEX(CELL_INDEX(I,J  ,K), 2)>0 .OR. WALL_INDEX(CELL_INDEX(I,J,K+1), 2)>0)  CYCLE
         CALL RANDOM_NUMBER(RN2)
         RN=REAL(RN2,EB)
         RN = NOISE_VELOCITY*(-1._EB + 2._EB*RN)*CELL_SIZE
         W(I,J,K)   = W(I,J,K)   - RN*RDY(J)
         W(I,J+1,K) = W(I,J+1,K) + RN*RDY(J+1)
         V(I,J,K)   = V(I,J,K)   + RN*RDZ(K)
         V(I,J,K+1) = V(I,J,K+1) - RN*RDZ(K+1)
      ENDDO
   ENDDO
ENDDO
DO K=1,KBM1
   DO J=1,JBAR
      DO I=1,IBM1
      !  IF (WALL_INDEX(CELL_INDEX(I,J  ,K), 3)>0 .OR. WALL_INDEX(CELL_INDEX(I+1,J,K), 3)>0 .OR. &
      !      WALL_INDEX(CELL_INDEX(I,J  ,K), 1)>0 .OR. WALL_INDEX(CELL_INDEX(I,J,K+1), 1)>0)  CYCLE
         CALL RANDOM_NUMBER(RN2)
         RN=REAL(RN2,EB)
         RN = NOISE_VELOCITY*(-1._EB + 2._EB*RN)*CELL_SIZE
         W(I,J,K)   = W(I,J,K)   - RN*RDX(I)*R(I)*RRN(I)
         W(I+1,J,K) = W(I+1,J,K) + RN*RDX(I+1)*R(I)*RRN(I+1)
         U(I,J,K)   = U(I,J,K)   + RN*RDZ(K)
         U(I,J,K+1) = U(I,J,K+1) - RN*RDZ(K+1)
      ENDDO
   ENDDO
ENDDO
DO K=1,KBAR
   DO J=1,JBM1
      DO I=1,IBM1
      !  IF (WALL_INDEX(CELL_INDEX(I,J  ,K), 2)>0 .OR. WALL_INDEX(CELL_INDEX(I+1,J,K), 2)>0 .OR. &
      !      WALL_INDEX(CELL_INDEX(I,J  ,K), 1)>0 .OR. WALL_INDEX(CELL_INDEX(I,J+1,K), 1)>0)  CYCLE
         CALL RANDOM_NUMBER(RN2)
         RN=REAL(RN2,EB)
         RN = NOISE_VELOCITY*(-1._EB + 2._EB*RN)*CELL_SIZE
         V(I,J,K)   = V(I,J,K)   - RN*RDX(I)
         V(I+1,J,K) = V(I+1,J,K) + RN*RDX(I+1)
         U(I,J,K)   = U(I,J,K)   + RN*RDY(J)
         U(I,J+1,K) = U(I,J+1,K) - RN*RDY(J+1)
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE INITIAL_NOISE


SUBROUTINE UVW_INIT(NM,FN_UVW)

! Read UVW file

USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER,SHUTDOWN
INTEGER  :: I,J,K,II,JJ,KK,IW,IOR,LU_UVW,IERROR,IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
INTEGER, INTENT(IN) :: NM
CHARACTER(80), INTENT(IN) :: FN_UVW
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
CHARACTER(3) :: STR
TYPE(WALL_TYPE), POINTER :: WC=>NULL()

CALL POINT_TO_MESH(NM)

LU_UVW = GET_FILE_NUMBER()
OPEN(UNIT=LU_UVW,FILE=FN_UVW,FORM='FORMATTED',STATUS='OLD',IOSTAT=IERROR)

IF (IERROR/=0) THEN
   WRITE(STR,'(I3)') NM
   MESSAGE = 'ERROR: Problem with MESH '//TRIM(ADJUSTL(STR))//'; UVW_FILE '//TRIM(FN_UVW)//' does not exist.'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
   RETURN
ENDIF

IF (PERIODIC_TEST==2) THEN
   IMIN = 1
   IMAX = IBAR
   JMIN = 1
   JMAX = JBAR
   KMIN = 1
   KMAX = KBAR
ELSE
   READ(LU_UVW,*) IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
   IMIN = MAX(0,IMIN)
   IMAX = MIN(IBAR,IMAX)
   JMIN = MAX(0,JMIN)
   JMAX = MIN(JBAR,JMAX)
   KMIN = MAX(0,KMIN)
   KMAX = MIN(KBAR,KMAX)
ENDIF
DO K=KMIN,KMAX
   DO J=JMIN,JMAX
      DO I=IMIN,IMAX
         READ(LU_UVW,*,IOSTAT=IERROR) U(I,J,K),V(I,J,K),W(I,J,K)
         IF (IERROR/=0) THEN
            U(I,J,K)=0._EB
            V(I,J,K)=0._EB
            W(I,J,K)=0._EB
         ENDIF
      ENDDO
   ENDDO
ENDDO

CLOSE(LU_UVW)

IF (PERIODIC_TEST==2) THEN
   U(0,:,:) = U(IBAR,:,:)
   V(:,0,:) = V(:,JBAR,:)
   W(:,:,0) = W(:,:,KBAR)
ENDIF

US=U
VS=V
WS=W

! Set normal velocity on external and internal boundaries (follows divg)

DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
   WC  => WALL(IW)
   IOR = WC%ONE_D%IOR
   II  = WC%ONE_D%II
   JJ  = WC%ONE_D%JJ
   KK  = WC%ONE_D%KK
   SELECT CASE(IOR)
      CASE( 1)
         WC%ONE_D%U_NORMAL_S = -U(II,JJ,KK)
      CASE(-1)
         WC%ONE_D%U_NORMAL_S =  U(II-1,JJ,KK)
      CASE( 2)
         WC%ONE_D%U_NORMAL_S = -V(II,JJ,KK)
      CASE(-2)
         WC%ONE_D%U_NORMAL_S =  V(II,JJ-1,KK)
      CASE( 3)
         WC%ONE_D%U_NORMAL_S = -W(II,JJ,KK)
      CASE(-3)
         WC%ONE_D%U_NORMAL_S =  W(II,JJ,KK-1)
   END SELECT
   WALL(IW)%ONE_D%U_NORMAL = WALL(IW)%ONE_D%U_NORMAL_S
ENDDO

END SUBROUTINE UVW_INIT

END MODULE INIT
